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Abstract. We introduce a class of convex, higher-dimensional billiard models which 
generalise stadium billiards. These models correspond to the free motion of a point-particle 
in a region bounded by cylinders cut by planes. They are motivated by models of particles 
interacting via a string-type mechanism, and confined by hard walls. The combination of these 
elements may give rise to a defocusing mechanism, similar to that in two dimensions, which 
allows large chaotic regions in phase space. The remaining part of phase space is associated 
with marginally stable behaviour. In fact periodic orbits in these systems generically come 
in continuous parametric families, associated with a pair of parabolic eigen-directions: the 
periodic orbits are unstable in the presence of a defocusing mechanism, but marginally stable 
otherwise. By performing the stability analysis of families of periodic orbits at a nonlinear 
level, we establish the conditions under which families are nonlinearly stable or unstable. As 
a result, we identify regions in the parameter space of the models which admit non-linearly 
stable oscillations in the form of whispering gallery modes. Where no families of periodic 
orbits are stable, the billiards are completely chaotic, i.e. the Lyapunov exponents of the 
billiard map are non-zero. 
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1. Introduction 

Billiard models, in which a point particle moves uniformly until it undergoes abrupt 
elastic collisions with a fixed boundary, are the playground of statistical physicists and 
mathematicians alike, whose work focuses on the interplay between their dynamical and 
statistical properties (T). 

There are two main categories of chaotic billiards. The better-known ones are (semi- 
)dispersing billiards, of which the hard-sphere gas is the prototypical example. The motion of 
hard spheres in a bounded region, undergoing elastic collisions with each other, is equivalent 
to that of a billiard model - a point particle which moves uniformly in the exterior of a 
collection of spherical cylinders in a high-dimensional phase space, with specular collisions at 
the boundary [0. The Sinai billiard [3] and the periodic Lorentz gas [4] are two-dimensional 
examples of dispersing billiards, in which a particle moves outside a disk on the 2-torus or a 
periodic configuration of them on the plane. The latter is a useful model of tagged particle 
diffusion in a binary mixture 0; it enjoys fast decay of correlations and thus converges to 
a Brownian motion [6j. The mechanism giving rise to chaos in these billiards, as well as 
in hard sphere gases, is that of dispersion, where nearby trajectories separate exponentially 
fast as a function of the number of collisions with surfaces of positive curvature, implying 
that the system has a positive Lyapunov exponent J7]] and is thus chaotic. Indeed, hard- 
sphere gases are essentially the only systems of interacting particles for which rigorous results, 
ranging from hyperbolicity to exponential decay of correlations, have been firmly established 

The second category is that of defocusing billiards, the paradigm of which is the 
Bunimovich stadium [QT1CE2]]. Here, chaos is due to a mechanism different from dispersion, 
namely the defocusing mechanism. As opposed to the Sinai billiard, the boundary of the 
stadium curves inwards with respect to the particle motion, i.e. it is convex. Although nearby 
trajectories initially focus after a collision with this boundary, if the distance to the next 
collision is longer than the distance to the focal point then they eventually defocus even 
more. This mechanism thus leads to an overall expansion in phase space, again measured 
by a positive Lyapunov exponent [ 13ll . 

Since its discovery, the mechanism of defocusing has attracted much attention in the 
physics community, particularly in connection to quantum chaos lfT4ll . acoustic experiments 
in closed chaotic cavities lfT5ll . optical microcavity laser experiments [[161 S3, or quantum 
conductance experiments lfT8ll to name but a few. 

In spite of its potential appeal to a broad range of physical applications, it has, 
however, remained a difficult problem to establish the conditions under which the defocusing 
mechanism can be extended to higher dimensions lfT9ll20ll . There are still only a few models 
of higher-dimensional stadia based on three- and higher-dimensional cavities which have lent 
themselves to a systematic study and are known to be fully chaotic B^l22l[23l[24ll25ll . 

There are arguably two difficulties which must be dealt with in order to obtain chaotic 
behaviour: (i) two curved dimensions, as in a spherical cap, often generate stable oscillations; 
and (ii) three-dimensional billiard domains with a single curved dimension, such as a 
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cylindrical surface, have flat components which complicate the stability analysis. 

The first of these two difficulties is related to the fact that one cannot simply construct a 
three-dimensional surface of revolution by rotating, for example, a two-dimensional stadium 
to obtain a chaotic three-dimensional cavity, since angular momentum is then conserved. 
Furthermore, as shown by Wojtkowski If20l . intersecting three-dimensional spherical caps and 
flat components may produce stable periodic orbits, giving rise to billiards with mixed phase 
space. The models studied by Bunimovich and Rehacek GH [22l [23) provide an exception to 
this principle and, in this sense, are rather special. They are furthermore non-convex. 

A natural alternative to spherical caps is to use boundary components made out of 
cylindrical caps and flat planes. It is clear that "extruding" or extending the two-dimensional 
stadium to a three-dimensional cylindrical-type stadium (a cylindrical shape whose cross 
section is a stadium) by a perpendicular translation does not suffice to produce a fully chaotic 
billiard in three dimensions, since the motion in this perpendicular direction is trivial. Thus 
Papenbrock [26] proposed to consider a three-dimensional stadium with half-cylindrical caps 
along perpendicular axes at both ends of a cuboidal shape. This construction gives rise 
to a billiard which is fully chaotic, as was recently shown rigorously by Bunimovich and 
Del Magno J24l [251, and is a particular case of more general constructions discussed by 
Wojtkowski |[27ll . In this respect, it is the first known example of a three-dimensional billiard 
in a convex domain with this property. 

The specificity of the Papenbrock three-dimensional stadium is that its cylindrical caps 
are placed opposite each other in such a way that periodic orbits visiting the caps are isolated. 
That is, even though each cylindrical component is curved along a single direction, , there is 
no degree of freedom with which to move the periodic orbits around, since the axes of the two 
cylinders are mutually transverse. 

The Papenbrock stadium thus avoids the second of the two difficulties we alluded 
to above, which are typical of three-dimensional cavities constructed out of cylindrical 
components and flat planes: namely, that their periodic orbits generically belong to continuous 
parametric families. Take, for example, the three-dimensional cylindrical-type stadium 
discussed above and break its translational symmetry by cutting it with an oblique plane. A 
periodic orbit of this system can be shifted along the direction of the cylinder axis, remaining 
basically unchanged, thus giving rise to a one-parameter family. This has the consequence 
that its stability analysis yields, out of the four phase-space dimensions of the billiard map, 
two parabolic directions (with unit eigenvalue), corresponding to motion along the family. 

Another example of a three-dimensional cavity whose periodic orbits belong to 
continuous parametric families is the three-dimensional billiard whose domain consists of 
the intersection of a sphere with a cuboid ll28ll . Though this billiard was numerically found to 
have four non-zero Lyapunov exponents in a large interval of its parameter, this property is yet 
unproven as it does not lend itself to the kind of analysis performed in refs B^[22ll23ll24l[25ll . 

It is the goal of this paper to consider billiards whose phase space combines neutral 
directions with curved regions, associated with stable and unstable behaviors, and explore the 
conditions under which they can display hyperbolic regimes, i.e. such that all pairs of opposite 
Lyapunov exponents are non-zero. To this end, we introduce a class of higher-dimensional 
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convex billiards, in which the particle motion occurs inside cylindrical domains bounded by 
oblique planes; we call these cylindrical stadium billiards. The reasons for using the term 
stadium are manifold. As we shall show, these billiards share some of the essential properties 
of the two-dimensional stadium billiard. In particular, and foremost, the mechanism that 
drives instability is of the defocusing type; other common features of specific interest to us, 
and which will be discussed below, are the existence of classes of periodic orbits similar to 
the bouncing-ball and whispering-gallery modes of the stadium billiard [291 . 

The motivation for such a general class of models arises from models of particles which 
interact via virtual strings, as introduced by Papenbrock to model self-bound nuclei J30l |3H 
and which were subsequently modified by Gilbert and Lefevere to model heat conduction 
by gas particles trapped in a nano-porous solid matrix ll32l . The interaction between pairs 
of neighbouring particles is a sort of hard-core interaction at a distance, and works as 
follows: particles move freely until they reach a critical relative distance, at which point 
they interchange the longitudinal components of their momenta. This can be thought of as 
an interaction mediated by a string, which affects the particles only once the string becomes 
fully extended, at which point it applies an instantaneous impulse to each particle so they 
flip directions. A chain of particles with such interactions were considered in [32], but placed 
inside a rectangular domain divided up into square cells, with each particle confined to its own 
cell. As their study shows, the combination of specular reflections along flat walls and circular 
arcs gives rise to chaotic motion in this system. Thus, given a system of N > 2 particles, they 
find 2N — 1 pairs of non-zero Lyapunov exponents opposite each other (the single pair of 
vanishing exponents corresponds to the energy conservation and associated time translation 
symmetry). 

As we shall shortly describe further, models of particles with such a string-type 
interaction, and confined by hard walls, are equivalent to billiards in higher-dimensional 
configuration spaces inside convex regions bounded by cylinders and planes. This motivates 
the introduction of a general class of convex billiards bounded by cylinders and planes, which 
includes such models. Note that this is the opposite of hard-sphere gases, which are equivalent 
to billiards outside cylinders. 

The presence of a flat direction along the axes of the cylinders in such systems 
implies that periodic orbits come in continuous parametric families, with each periodic orbit 
displaying marginal stability along the direction of the family (each dimension associated with 
a flat direction actually gives rise to a pair of parabolic eigenvalues along the corresponding 
direction and conjugate momentum). Although the existence of such a family was noted in 
ref. [33], its detailed analysis is notably absent in the models previously considered. Indeed, 
the stability of such families of periodic orbits must be studied at a nonlinear level in order to 
determine whether or not part of the family can or cannot be stabilised - a daunting task. 

As a first step towards a comprehensive understanding of the dynamics of such systems, 
we present in this paper a detailed analysis of the dynamical regimes occurring in one of the 
simplest such systems, which consists of a cylinder in three dimensions with a flat bottom and 
cut by an oblique plane at its top. This analysis will be performed in terms of the stability of 
periodic orbits, and comes in two parts. At a linear level, stability is understood in terms of 
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the lengths of the segments of the orbit, and, specifically, that of paths corresponding to two 
collisions with the cylindrical surface, separated by a single collision with the oblique plane. 
As we shall show, when this length is large enough, defocusing takes place and the periodic 
orbit is unstable. If on the contrary, this length remains small, then the orbit may remain 
stable. However, because such periodic orbits have a pair of parabolic eigenvalues, the effect 
of nonlinear perturbations along the neutral directions must be taken into account, which is 
the second part of the stability analysis. 

There are essentially two classes of periodic orbits of cylindrical stadium billiards whose 
segments' lengths can remain small enough that they are potentially stable. We identify these 
two classes as planar and helical periodic orbits. The former are similar to bouncing-ball 
orbits in the stadium billiard, and the latter to whispering-gallery modes. By analysing in 
detail the existence conditions and linear and nonlinear stability properties of these two classes 
of orbits, we show, at the linear level, the existence of a bifurcation from marginally- stable to 
unstable (hyperbolic) regimes, and, at a nonlinear level, the existence of a "restoring force" 
which, in specific regions of the parameter space, is able to stabilise the motion by restricting it 
to a certain part of the parameter space where the family of periodic orbits remains marginally 
stable, thus preventing the periodic orbit from escaping to a hyperbolic region through the 
bifurcation point. The parameter regions where none of the periodic orbit families can be 
stabilised correspond to chaotic regimes of the system. 

The paper is organised as follows. In section [2} we provide a detailed description 
of particles interacting through a string-type mechanism in terms of cylindrical billiard 
models. In section[3]we introduce the particular three-dimensional cylindrical stadium billiard 
described above, and we offer a general characterisation of its periodic orbits in section]?} The 
two classes of planar and helical periodic orbits are studied in further details in sections [5] and 
[6j respectively, where we identify parametric regions of non-linear stability of these periodic 
orbits. A summary and discussion of our results are presented in section [7} Two appendices, 
Appendix~~A| and |Appendix~l3} provide exact forms obtained for the linear and non-linear 



analysis of some of the periodic orbits studied in sections [5] and [6j 



2. Interacting-particle models and cylindrical stadium billiards 

We start by relating interacting-particle models with string-type interactions and billiard 
models, and then introduce the general class of cylindrical stadium billiards. 



2.7. Particles with string-type interactions and flat walls 

Consider N >2 particles interacting via a string-type interaction [[30l[3ll[32l. The particles 
are joined by strings in some configuration, and have the following dynamics. Two particles 
which are joined by a string move freely until they reach a mutual distance r/y = £ 9 at which 
point the string becomes taut (fully extended), and exerts an instantaneous impulse on each 
particle, which interchanges the components of their momenta along the line between them. 
The interaction potential between the particles due to the string, as a function of the distance 
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r, is thus v(r) = for r < £ and v(r) = °° for r > £ 

The dynamics of the A/-particle system corresponds to the motion of a point F = (q, p) = 
(qi , . . . , q#, Pi , . . . , Pw) in a 2dA/-dimensional phase space, where d is the spatial dimension 
of the system, so that q/ G R d and p/ G K J . 

If particles / and 7 are joined by a string of length £, then they are constrained by the 
inequality 

l|q/-q;l| 2 <^ (D 

which specifies a hyper- spherical region S G R 2d in the dimensions spanned by the coordinates 
qi and q 7 . The other coordinates are free, so that this corresponds to an allowed cylindrical 
region S x R d ( N ~ 2 ) in configuration space. This is analogous to the case of hard-sphere 
interactions J2|, for which the constraint has the opposite inequality, ||q/ — q 7 -|| 2 > £ 2 . 




Figure 1. Two interacting particles joined by a string and trapped in a square box. Their 
relative positions are constrained to a disk of radius £, while each particle is reflected elastically 
by the horizontal and vertical walls. The two continuous lines show the trajectories of each 
particle as they interact and collide with the walls. The motion of the center of mass (dashed 
line) is uniform in between collisions of either particle with the walls. The dotted arrows 
indicate the directions of the velocity vectors at the initial and final positions. 

The particles may also be confined by hard walls such as shown in figure [TJ where two 
particles joined by a string are trapped to a square box. Each wall imposes a further constraint 
on the position q; of each particle, namely (q; — c) • n < 0, where c is a point on the wall and n 
the outgoing normal vector. Each wall thus imposes N constraints on the total configuration- 
space position vector q, which are of the form (q — c) • n < 0, where A = (0, . . . , n, . . . , 0), with 
n appearing in the ith position, and similarly for c; these correspond to half-spaces bounded 
by hyper-planes. 

The configuration space vector q is thus confined to the intersection of the cylinders 
coming from the string- type interactions, and the half- spaces corresponding to the bounding 
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hard walls; this region is the intersection of convex sets, and is thus itself convex. Due to 
the fact that the interactions are hard-core, the motion in phase space is in fact equivalent 
to a semi-focusing billiard inside this convex set [24]; this should be compared to the semi- 
dispersing billiards which arise from hard- sphere gases 0. 

2.2. Simplest interacting -particle model 

In order to get a grasp on the types of phenomena which can occur in such models, it is 
instructive to look at the simplest model in this class, which consists of a two-dimensional 
string-bound diatomic molecule, with string length (maximum separation) £, confined in 
a channel formed by two parallel, infinitely long, hard walls at y = ±|, i.e. similar to 
the example shown in figure [TJ but with the vertical walls removed and horizontal walls 
infinitely long. We denote the particle positions and velocities by q/ = (*/,)>/) and p/ = (wj, v/), 
respectively, for i — 1,2. The string constrains the particles by ||qi — 1 1 2 < and the 
above argument shows that the hard walls give 4 allowed half-spaces, given by extending the 
inequalities — f < J; < f to the configuration space vector q. 

Although the phase space of this system is a priori 8-dimensional, translational symmetry 
along the channel implies that the horizontal component of linear momentum of the centre of 
mass is conserved, so that we may choose a reference frame in which it vanishes, u\ + U2 = 0, 
and in which we can also fix its position, x\ + x^ — 0. It is convenient to rotate the 
coordinate system by angles n/4 and normalise them according to £, setting X = (x\ — X2)/£ 9 
Y = (yi — yi)/^, and Z = {y\ -\-y2)j£. This gives an equivalent billiard system, consisting of a 
point particle which bounces elastically inside the region bounded by the cylinder X 2 + Y 2 < 1 , 
and chopped by the four planes Y ±Z = ±h 9 with h = djl. This, in turn, is equivalent to a 
billiard in a domain bounded below by the plane Z = and above by the two planes Z± Y = h, 
shown in figure [2j 

Similar constructions give rise to cylindrical stadia in higher dimensions. For example, 
in the case of the diatomic molecule inside a 2-dimensional square box with hard walls 
shown in figure [TJ there is no longer conservation of momentum, so that we obtain a billiard 
inside a cylinder with two curved and two flat directions in 4 dimensions, bounded by eight 
hyperplanes. 

2.3. Cylindrical stadium billiards 

The examples of particles interacting via string-type interactions motivate our introduction 
of a large class of cylindrical stadium billiards. Informally, these consist of one or more 
cylinders which are cut by (hyper-)planes. More formally, they are the intersection of one 
or more infinite cylindrical region(s) with two or more half spaces, such that the resulting 
billiard is bounded and such that at least part of the boundary of the resulting region is curved 
(that is, comes from a cylinder). This last condition is required since otherwise the region is 
polyhedral, and the dynamics is known to be non-chaotic lfT3ll . 

Note that these billiards do not share the property of the Bunimovich and Papenbrock 
stadia that the cylindrical and flat portions of the boundary meet tangentially (C 1 boundaries). 
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Figure 2. The motion of two interacting particles joined by a string in a two-dimensional 
infinite channel is equivalent to the motion of a point particle in a three-dimensional cylindrical 
domain chopped by four planes, or, equivalently, a cylindrical domain with a flat bottom and 
topped by two interesting planes at z =b y = h, as shown above. This billiard has a single 
parameter, the height of the oblique planes, given by the ratio of the width of the channel to 
the length of the string of the interacting di-atomic model. 

Nonetheless, they are higher-dimensional convex billiards with large chaotic regions, and in 
this sense are a possible generalisation of the stadium which are of intrinsic interest. Note 
also that the stadium is only one of a large class of 2D billiards proved to be chaotic, many of 
which do not share its property of having C 1 boundaries [T3j. 

3. A three-dimensional cylindrical stadium billiard 

In order to investigate in more detail the types of mechanisms which are present and which 
give rise to chaotic and/or stable motions in cylindrical stadium billiards, we henceforth 
restrict attention to one among the simplest possible such systems, which is a simplified 
version of the 3D cylinder found in the previous section. 

We consider a cylindrical stadium billiard in M 3 formed by the cylinder x 2 + y 2 < 1 with 
unit radius and vertical axis, and bounded below by the plane z = 0, which is perpendicular 
to the axis, and an oblique plane that is inclined away from perpendicular. Here we focus 
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on a one-parameter family of such billiards, where the inclined plane is placed at an angle 
k/4 with respect to the cylinder axis, and is parameterised by its height h, with equation 
y + z = 1 + h. The parameter h> —2 thus defined measures the minimal height of the cylinder 
at y = 1 ("above" its base); h < corresponds to the case where the inclined plane cuts the 
flat plane inside the cylinder. The convex billiard domain then consists of the intersection of 
the cylinder with the two half-spaces z > and y + z< l+h-stt figure [3} 

We view this convex billiard as the elementary cell of an expanded square- shaped 
cylindrical billiard, or more simply squared cylindrical stadium, obtained by "unfolding" the 
elementary cell, that is, repeatedly reflecting it in its two planes. This cavity is thus made 
out of the union of four sections of cylinders, which intersect pairwise at right angles - see 
figure |4| The only relevant parameter is the height of each of these cylinders; their half-height 
is the parameter h introduced above, which we will refer to as the geometrical parameter. It 
is bounded below by h — —2, but is not bounded above. 




Figure 3. Examples of the cylindrical billiard domains we study: (a) h = — 1; (b) h = 0; (c) 
h=\. 




Figure 4. Squared cylindrical stadia obtained by repeatedly reflecting the cylindrical billiard 
domains displayed in figure|3]in their bottom and top planes, (a) h = — 1; (b) h — 0; (c) h = 1. 
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Considering the unfolded billiard, we will refer to the plane generated by the axes of the 
four cylinders as the plane of the billiard. A transverse plane refers to any plane perpendicular 
to that plane. These include the planes perpendicular to the cylinders' axes, as well as the 
oblique planes in which they intersect. 

Choosing a Cartesian coordinate system whose center is at the center of symmetry of 
the square, we take the plane x = to be the plane of the billiard. The cylinder axes are then 
located at x = and y = ±(l+h), parallel to the z-axis, and at x = and z = ± ( 1 + h) , parallel 
to the y-axis. 

The unfolded billiard undergoes significant topological changes at h — — 1 and h — 0. 
When h = — 1, each pair of parallel axes coalesces, so that the unfolded billiard is a convex 
cavity for — 2 < h < — 1. When = 0, two parallel cylinders intersect along a single line on 
their surfaces - they intersect for h < and separate for h > 1 . The unfolded billiard thus has 
genus for h < and genus 1 for h > 0. 

4. Periodic orbit families 

Much in the same way that one distinguishes bouncing-ball and whispering-gallery periodic 
orbits in the two-dimensional stadium billiard, we will distinguish two specific types of 
periodic orbits in the three-dimensional square cylindrical stadium. The first type consists 
of planar periodic orbits, which are analogous to the bouncing-ball periodic orbits: they 
bounce from one side of the cavity to the opposite one, within the plane of the billiard. The 
second type consists of helical orbits, which whirl around the surfaces of the cylinders, going 
smoothly from one cylinder to the next across their intersection. As we shall show, the analysis 
of these two classes of periodic orbits is sufficient to determine the parametric regions where 
stable oscillations in cylindrical stadium billiards occur. 

The specificity of cylindrical billiards is that all their surfaces have a flat component. 
This has the consequence that every periodic orbit has at least two parabolic eigen-directions, 
i.e. whose corresponding eigenvalues are unity. These are associated to motion parallel to the 
axes of the cylinders. 

Since the squared cylindrical stadium gives rise to a four-dimensional symplectic billiard 
map, a periodic orbit also has another pair of eigenvalues, in addition to these two parabolic 
ones, whose product is also unity. This pair is associated to the motion in planes transverse 
to the cylinders' axes. It can be either hyperbolic, when the two eigenvalues are real and are 
inverses of each other, or elliptic, when they are complex conjugates with unit modulus. The 
case with four parabolic eigenvalues is irrelevant. 

The consequences of the existence of a pair of parabolic eigenvalues for each periodic 
orbit are twofold. First, a periodic orbit is unstable when the second pair of eigenvalues 
is hyperbolic, but only marginally stable when this pair is elliptic. In other words, a 
given periodic orbit can never be linearly stable by itself. The second consequence is that 
periodic orbits exist in continuous one-parameter families, which can be parameterised by 
a displacement along the cylinders' axes. We refer to this displacement as the dynamical 
parameter, which we will denote in general by £, and which is defined suitably to parameterise 
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each family of orbits. It is not to be confused with the height of the cylinders, which is the 
model's geometrical parameter h. 

If one were to limit the stability analysis to linear terms, as has previously been the 
case, then the naive conclusion would be that a family of orbits is stable on condition that 
every single orbit of this family admits a pair of elliptic eigenvalues, and provided that there 
are no geometric constraints which prevent this family from being complete because of a 
discontinuity. The reasoning is that the parabolic eigen-directions behave neutrally, so that 
given a periodic orbit and a small perturbation along these directions, the trajectory will move 
along them as it oscillates in the plane spanned by the elliptic eigenvectors, sweeping through 
the whole family of associated periodic orbits, parameterised by different values of £, but 
left unscathed in the absence of a destabilizing mechanism. If, on the other hand, the family 
of periodic orbits undergoes a transition from elliptic to hyperbolic regimes as £ varies - 
as happens in most cases - then the initially stable oscillations in the elliptic regime would 
become unstable at the bifurcation, when a value of £ is reached at which the family crosses 
over from the elliptic to the hyperbolic regions, and thus the motion would destabilise. 

As it turns out, when elliptic regions occur for a given family of periodic orbits, 
transitions from elliptic to hyperbolic regimes are always found as soon as the value of the 
maximal height 2 + h of the cylinders in the elementary cell is larger than the cylinders' radii, 
i.e. when h > — 1. The reason is essentially that beyond a bifurcation value of the dynamical 
parameter associated to a given family of periodic orbits, the distance between successive 
collisions is, at some point on the orbit, large enough to induce defocusing. According to the 
above arguments, one would then be led to believe that no stable oscillations could persist for 
square cylindrical stadia with heights larger than this bound. 

However, interestingly, the story does not stop there. Extending the stability analysis 
to the next order beyond linear order, we find that, under some assumptions on the initial 
conditions, nonlinearities are sometimes able to stabilise a (partial) family of periodic orbits in 
a limited region of the dynamical parameter £, leaving the orbit confined to the elliptic region, 
oscillating back and forth between two values of this parameter. Remarkably, the nonlinear 
perturbations are quadratic in the coordinates associated to the elliptic eigenvalues. As we 
show below, the key to understanding whether there exists a nonlinear stabilizing mechanism 
lies in a qualitative analysis of the eigenvalue spectrum of these quadratic forms: the signs 
of their eigenvalues tell us whether or not the oscillations in the plane spanned by the elliptic 
eigenvectors are able to stabilise the oscillations along the parabolic directions. 

As our analysis will demonstrate, the conditions for stable oscillations are met only for 
a restricted set of values of the geometrical parameter h, consisting of a disjoint union of 
intervals and isolated values. For values of the parameter in that set, the phase space is 
typically mixed - different stable regimes may coexist, together with unstable ones. On the 
other hand, away from these values, the square cylindrical stadia are fully chaotic. 



Higher-dimensional convex billiards with cylindrical shape 



12 



5. Planar orbits 

By planar orbits, we mean those that remain in the plane of the billiard (x = 0), corresponding 
to all the orbits whose initial conditions are in that plane, with velocity components along the 
y- and z-axes onljjf] 

For h < 0, the periodic orbits in that plane are readily identified as the periodic orbits 
of a square billiard with sides of length 2(2 + h). Given the relative prime integers m and n, 
the point q = {0 5> y, z} is periodic provided the velocity is given by p = {0, cos 0) m ^, sin 0) m ^}, 
where (O m ^ n = arctan(fz/m). Over a period, the corresponding orbit makes exactly m collisions 
on each of the horizontal walls, and n on each of the vertical walls, in a time x m ^ n — 
4(2 + h)Vm 2 + n 2 max(m/n,n/m). 

Furthermore, for any pair of relative primes {m, n}, there is a continuous family of orbits, 
indexed by these integers, obtained by varying the position of the collision impacts on the 
square sides. Still assuming h < 0, this family extends over all the possible values of y and z, 
viz. -(2 + h) <y,z<2 + h. 

For h > 0, the periodic planar orbits are those of a two-dimensional square billiard of 
side length 2(2 + h) with a square of side length h removed in the centre. The periodic orbits 
of this billiard include some that do not collide with the central square, which are common 
to the h < case, as well as others that do collide with the central square and which are in 
general more difficult to classify. 

We will call (m,n)o planar periodic orbits (PPO) those periodic orbits of the plane that 
make collisions only with the outer surfaces of the cylinders. Examples are shown in figure [5] 
Except for n = and m — n — 1 , we will assume m <n. 
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Figure 5. The periodic orbits on a square domaine are specified by the direction of the velocity 
vectors with slopes tan(^/m), where m : n G N. Several orbits of the (1, l)o and (l,2)o families 
are shown, corresponding to the height parameter h = 3/4. 

$ There are other planar orbits in planes transverse to thecylinders axes, such as a period-four orbit consisting 
of a bowtie figurebetween two opposite circlular arcs and passing through the center of symmetry of thebilliard. 
Such orbits exist for — 1 < h < and can be found near the y = and z = planes, but are not significant since 
they are always unstable. 
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5. 7. ( 1 , 0)o planar periodic orbits 

5.1.1. Existence The simplest class of planar orbits are period-2 orbits which bounce back 
and forth between two opposite cylinders with velocity perpendicular to the cylinder axes. 
Assuming they propagate along the y axis, we identify them by the phase point with Cartesian 
coordinates 



(2) 



These orbits exist everywhere along the axis of the cylinder so long as h < 0, i.e. 
— 2 — h < £ < 2 + h, or away from the center when h > 0, i.e. h < \e\ <2 + h. Specific 
examples are displayed in figure [6} 





( x 


\ 




( 






/ 


UQ 


\ 




( 





\ 


qo = 


yo 








-2-h 


•> Po = 




vo 








1 






\ zo 


) 




\ 


e ) 




V 


Wq 


J 




{ 





J 



h=-l,5 




(a) 




(b) 



(c) 



Figure 6. Period-2 (l,0)o PPOs of the square cylindrical stadium, (a) h = — 3/2, (b) h = — 1, 
(c) h = — 1/2. These orbits correspond to £ = in but similar orbits exist of all possible 
values of £ along the vertical axis. 



5.7.2. Linear stability The (l,0)o PPOs are stable only in the interval —2<h< —1. This 
result is a straightforward transposition of the stability of periodic orbits bouncing back and 
forth between two opposite circular arcs of identical curvatures: h = — 1 corresponds to the 
case where the distance between the arcs is equal to their diameter. Thus the case (a) in 
figure [6] is stable, case (b) marginal and case (c) unstable. 

Cartesian coordinates. It is useful to go through the stability analysis explicitly. To do so, 
we compute the Jacobian of the periodic orbit by following the transformation rules of the 
tangent vectors 5q f and <5p; along the orbit, where / = x,y,z and each of these vectors has six 
dimensions. 

These transformation rules are given according to the two phases of propagation and 
collision by lf34ll 

H- * H + T5p ; , (propagation T)(3) 
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and 



(collision) (4) 



<5p, i — y Spi — 2 ySpi • nj n — 2 • Sn^ n — 2 ^p • <5n/ 

where n is the unit vector normal to the surface, and <5n is a triplet of six component vectors 
which depends on the curvature of the surface, 



p n 



(5) 



Given a cylinder with axis in the y- or z-directions and a trajectory in the plane of the billiard 
(the x components of the positions and velocities vanish), this reduces to 

8n x = 8^ , (6) 
8n y = 8n z = 0. (7) 

The orbit starting at the coordinates Q bounces back and forth from qo, po to 

qi = (0,2 + M), Pi = (0,-1,0). (8) 

Correspondingly, letting T = 2(2 + h) denote the time between successive collisions, 
the tangent vectors are mapped according to equations (|3j)— (|4j) as follows, following the 
propagation (prop.) and collision (coll.) phases: 
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coll. 
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1 -2t 








2t(1 -t) 
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4(1 -T) 








(1-2t) 2 -2t 
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(9) 



/ 



This matrix has six eigenvalues. Two pairs are trivial, the first corresponding to the 
conservation of energy and the associated time-translation symmetry, and the second being the 
pair of parabolic eigenvalues which are due to the flat components of the cylinders' surfaces. 
The only pair of non-trivial eigenvalues of this matrix are then 

Ai, 2 = 8/i(3 + h) + 17 ± (3 + 2h) y/(l+h)(2 + h). (10) 

They are elliptic for — 2 < h < — 1 and hyperbolic otherwise, as expected. 

We note that a similar conclusion (with simpler eigenvalues) can be reached by 
considering only half of the orbit. This is because the corresponding periodic orbit in the 
elementary cell involves only a single collision with the cylindrical surface. However, because 
the parity of the number of reflections on the transverse planes can be odd in the elementary 
cell, considering the periodic orbits in the expanded billiard is better suited to the non-linear 
stability analysis discussed below. 



Birkhoff coordinates. The symplectic structure of the collision map is best recovered by 
considering, instead of Cartesian coordinates, the Birkhoff coordinates of the billiard map, 
i.e. the coordinates on the constant energy surface in which the map is volume-preserving. 
They are given by: (i) 0, the position angle, measured in the plane transverse to the cylinders; 
(ii) £ , the sine of the associated velocity, measured with respect to the normal to the cylinder's 
surface; and (iii)-(iv) the pair z, w of the position and velocity measured along the cylinder 
axis: 

yo + 1 + h 71 

~~r 



arctan : 



x 



(lla) 



£o = sin ( arctan 

zo = e; 
0. 



= 0; 



w 



(lib) 
(He) 

(lid) 



In terms of the Birkhoff coordinates, the linear map ^ for the perturbations 80, 8%, 8w and 
8z reduces to 



f 80\ 

Sw 

V 5z I 



/ 




m et; 
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m^e 
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Sz J 



(12) 



It turns out that the structure of this linear map is common to all planar periodic orbits; the only 
non-trivial coefficients are the four matrix elements corresponding to the motion in the 0-^ 
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plane and the element of axial velocity contribution to the displacement along the cylinder 
axis. In the case of the (l,0)o PPO, the matrix of elements m aj g is given by 

/ ll + Sh(3 + h) -4(2 + h)(3 + 2h) 

4(l+h)(3 + 2h) ll + Sh(3 + h) 

1 

4(2 + /*) 1 J 

The four eigenvalues of this matrix form the pair of parabolic eigenvalues and the pair Ai 2 
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me 


m 6t; 





o\ 




m^ e 
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ITlzw 





o\ 






V 



(13) 



obtained in equation ( |T0| ). 

This result is summarised in figure [7| where the parameter space of the (l,0)o PPOs is 



colored, when the orbit exists, in white where the eigenvalues ( |T0| ) are elliptic and gray where 
they are hyperbolic, and in light red where the orbits do not exist. Similar conventions will be 
used throughout the paper. 




Figure 7. Results of the linear stability analysis of the (l,0)o PPOs plotted in the parameter 
space, with the dynamical parameter £ on the horizontal axes, and geometrical parameter h 
on the vertical axis. The light red areas represent the regions outside the range of allowed 
parameter values. The grey areas correspond to regions where the eigenvalues ( fT0| ) are 
hyperbolic, and white where they are elliptic. The line which appears in the elliptic region 
is a line of parabolic eigenvalues. 



The conclusion we infer from this calculation is that the parameter region — 2 < h < — 1 
corresponds to a regime with mixed phase space, where the period-2 orbits in the plane of the 
billiard give rise to an elliptic island. Though each of these orbits is only marginally stable, 
it belongs to a continuous family of similar orbits, which is moreover complete, meaning that 
it extends from one corner of the billiard to the opposite one. The family of orbits therefore 
corresponds to a plane of elliptic stability. 

We note that the (l,0)o PPOs are peculiar in the sense that their nonlinear stability 
analysis yields non-trivial corrections only to cubic order in the perturbations 86, 8^ and 
8w. Turning to the class of (l,n)o PPOs, we will see that, not surprisingly, the (l,0)o PPOs 
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are not the only planar periodic orbits which give rise to a stable regime. In fact, a perturbation 
of the (l,0)o PPOs in the plane of the billiard will typically be closer to a family of (l,n)o 
PPOs with n large. Interestingly, such periodic orbits may exhibit a mechanism for non-linear 
stability quadratic in the perturbations 89, 8% 9 <5w, absent for the (1,0) family. 

5.2. (1, l)o- planar periodic orbits 

5.2.1. Existence The orbit with unit speed at angle 71/ 4 with respect to the cylinder axes is 
the period-4 orbit which cycles through the phase points 

q = (0, -2 - h, e) qi = (0, -8, 2 + h) q 2 = (0,2 + h, -e) 
Po = ^ (0,1,1) ^ pi = ^(0,1,-1) ^ p 2 = -L (0,-1,-1) 



q 3 = (0,8,-2-/*) q 4 = qo 
1 

V2 

where -(2 + /i)<£<2 + /iif-2</i<0and/z-2<£<2-/iif0</i<2 



^ n WO 1 1 ^ * > (14) 

P3 = ^(0,-1,1) P4 = P0 



5.2.2. Linear stability Going through the analysis described in section [5X2} we obtain a 



linear mapping describing the evolution of the perturbation vector along the periodic orbit, 
which, in Birkhoff coordinates, is identical to equation ( [T2| ) with matrix elements 

m ee = 1 -4£ + S[h(2 + h) - £ 2 ][(1 +h) 2 - £ - £ 2 ], (15a) 

m e £ = - 4[2h(h + 2) -2e 2 + l][h(h + 3) - £ 2 + 2], (15b) 

m^ e = -4[2h(h + 2)-2£ 2 + l][h 2 + h-£ 2 ] (15c) 

m§§ =8[h(h + 2)-£ 2 ][(h+l) 2 -£ 2 + £]+4£ + l, (15d) 

= %V2(2 + h). (15c) 



m 



zw 



The two non-trivial eigenvalues are 

l l2 = 1 + 8 \-2h{h + 2)£ 2 + h{h + \f(h + 2) + £ 4 - £ l 



±4^(h-£ + l)(h + £ + l)[2h(h + 2)-2£ 2 + l] 2 [h(h + 2)-£ 2 ]. (16) 

The regions of elliptic and hyperbolic eigenvalues are displayed in figure [8} There are two 
separate elliptic regions, bounded by h = — 1 ± £ and h = — 1 ± y/l + £ 2 . For h = — 1 and 
2/3 < h < 2, there are no elliptic regions along the £ axis and, where they exist, the (1, l)o 
PPOs are everywhere hyperbolic for these values of the geometric parameter. 

The transitions from elliptic to hyperbolic regimes are understood as follows. The 
segment of orbit between, say, qi and q2 lies in the transverse plane z — y — 2 + h — £, which 
cuts the cylinders along the y- and z-axes in the form of two partial ellipses, glued one against 
the other. The semi-minor axis of these ellipses is unity (the radius of the cylinders) and the 
semi-major axis \[2. The segment of orbit that joins qi to q2 goes along the major axes of 
these two ellipses. The orbit is thus stable only so long as the distance traveled along the 
major axes is less than the radius y/l/A of the disk inscribed into the ellipses at qi and q2. 
This condition yields the critical values £ = ±(1 +h) 9 which separates the hyperbolic and 
elliptic regions in figure [8j 
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*s 




Figure 8. Results of the linear stability analysis of the (1, l)o PPOs plotted in the parameter 
space. The conventions are identical to figure [7] 



Our analysis shows that the family of (1, l)o PPOs under consideration is linearly stable 
and complete only for —2<h< —3/2. For larger values of h, though some orbits display 
elliptic behavior, they do not correspond to linearly stable oscillations in the vicinity of the 
periodic orbit. The reason, according to the linear stability analysis, is that a trajectory initially 
close to an elliptic PPO will oscillate in the planes transverse to the cylinder's axes, and 
simultaneously move along the z-axis due to perturbations along the w-axis, which remain 
unchanged under iterations of the map ( [T2| ), whose coefficient m zw , given by equation ( |15e| )> 
depends solely on h. In other words, the orbit moves in the parameter space of the (1, l)o 
PPOs along the £-axis due to the action of the perturbation along the vertical velocity axis. In 
the absence of nonlinear effects acting on Sw 9 the trajectory would eventually have to cross 
over from the elliptic to the hyperbolic region and therefore lose stability as the parameter £ 
translates away from the interval of elliptic eigenvalues. 



5.2.3. Nonlinear stability Linear stability would have us conclude that periodic orbit 
families whose non-trivial eigenvalues display bifurcations from elliptic to hyperbolic regimes 
must correspond to unstable oscillations, even when the reference periodic orbit we initially 
perturb away from is marginally stable. This, however, is in general incorrect. Going beyond 
linear stability, we find that the perturbations along the w and z coordinates are nonlinear 
functions of the oscillations that take place in the 6-% plane. At second order in the stability 
analysis, the 6-^ oscillations give rise to quadratic terms which act on the perturbations along 
the parabolic directions according to 

Sw^Sw + Q w (S0,S%), (17a) 
Sz ^ 8z + m zw 8w + Q z (8G,8^), (lib) 
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where Q W (S6, 8%) and Q Z (S6, 8%) are quadratic forms in their variables, which we write as 
Q W {89M)= (SO 8S) ( T: ^ )( « V (18fl) 




&(*M$)=(** 5§)( £ ^ | ( *® | . (18&) 
The coefficients <z w , Z^, c w , Z? z and c z are functions of the two parameters h and £. 




Now, if <5w is large enough to start with, then the quadratic terms in ( |17a| )-( [T7^| ) will 
do little to affect the motion along the w and z axes. Thus the problem we need to address is 
the following: assuming!] that 8w is small with respect to quadratic terms in 86 and 8%, we 
must find out whether the quadratic forms Q w and Q z are or are not able to induce oscillations 
of 8w and <5z, i.e. £. Such oscillations would prevent the contribution m zw 8w in (12) from 
inducing a translation in the parameter £ (along the z axis) and would thus keep the trajectory 
oscillating within the elliptic region. In this case the corresponding family of (1, l)o PPOs 
would remain stable as a result of the nonlinear oscillations. 

Fortunately, the solution to this problem is simpler than it may seem, for, since the 
coefficient m zw in ( |15e| ) is positive, a mechanism inducing oscillations of w around £ = 
(or, for that matter, around any fixed value £ = £o) exists whenever Q w has the opposite sign 
from £ (or £ — £q if £q ^ 0), i.e. whenever Q w is positive where £ is negative and Q w is negative 
where £ is positive. These oscillations in £ will be further amplified by Q z , constructively so 
provided Q z also takes the opposite sign from £. If, on the other hand, Q w is positive where 
£ is positive and is negative where £ is negative, then, whether or not Q z keeps the opposite 
sign from £, 8w will keep increasing or decreasing, steadily if Q w has the sign of £, or on 
average if it has the opposite sign from £, until the trajectory reaches the limits of the elliptic 
region and thus destabilises, either because of a transition to a hyperbolic regime, or because 
the limits of the range of allowed £ values is reached. 

The quadratic forms ( |18a| )-( [T8^| ) are positive- or negative-definite whenever the 
eigenvalues of their defining matrices have the same signs. (Since these matrices are 
symmetric and have real entries, their eigenvalues are always real.) Thus, irrespective of the 
oscillatory motion of the 0-t, variables, Q w (resp. Q z ) is positive when its two eigenvalues are 
positive and negative when they are negative. When the two eigenvalues are of opposite signs, 
on the other hand, we must bear in mind that the oscillations of the 6-^ variables keep going 
at their own pace, and they will probe the two eigendirections of Q w as they do so. Thus, as 
Q w (resp. Q z ) oscillates, it will be expected to take, on average, the sign of the eigenvalue that 
has the largest magnitude. 

Figure [9] shows the eigenvalue spectrum of the matrices of the quadratic forms ( |18a| )- 
( |18£| ), superimposed on the elliptic regions found from the linear stability analysis in figure [8| 
The coefficients of these quadratic forms are given explicitly in |Appendix A.l[ equations 
( |l.la| )-( [LTc| ) and ( |1.2a| )-( [L2c| ). For convenience, we analyze them by using a representation 
in four different colors, in terms of the sign and magnitude of their eigenvalues, which are 

§ Identifying the initial position along the z axis with the dynamical parameter £, we can always set Sz = 0. 



Higher-dimensional convex billiards with cylindrical shape 



20 




(a) -2 -1 1 2 (b) -2 -1 1 2 

Figure 9. Results of the nonlinear stability analysis of the ( 1 , 1 )o PPOs plotted in the parameter 
space. The conventions are similar to figure [8] with the elliptic regions filled according to 
the signs of the determinant and trace of the quadratic forms (a) Q w and (b) Q z \ in blue 
where both eigenvalues are positive (X1X2 > 0, X\ + X2 > 0), red where both are negative 
(Ai A-2 > 0, X\ + A2 < 0), green where one is more positive than the other is negative (X\ X2 < 0, 
X\ + A2 > 0), and cyan where one is more negative than the other is positive (X\Xi < 0, 
Ai+A 2 <0). 

denoted in both cases by X\ and X2. We observe that the elliptic regions below and above the 
lines at h = — 1 in both cases display clearly different patterns. Whereas the upper elliptic 
regions are everywhere unstable under the quadratic perturbations, the lower elliptic regions 
are stable around £ = 0. 




(a) (b) » (c) 

Figure 10. Nonlinearly stable oscillations of the (l,l)o PPOs measured at h = — 1.2. (a) 
Actual trajectory; (b) Linear oscillations measured in the 6-t; plane; (c) Nonlinear oscillations 
measured in the w-z plane (here we identify the z coordinate with the dynamical parameter £). 
Notice the order of magnitude of the £ oscillations as opposed to those of Sw. 



The separation between nonlinearly stable and unstable elliptic regions at h = — 1 is 
supported by numerical results such as those shown in figures [TO] and [TTj 
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Figure 11. Nonlinearly unstable oscillations of the (1, l)o PPOs measured at h = —0.1. (a) 
The iterations of the actual trajectory remain stable over a short time only, (b) The motion 
along the z axis presents a drift as the dynamical parameter £ decreases to the negative values 
and away from the elliptic regime, (c) The oscillations of 6 (as well as not shown) are 
initially elliptic, but soon hit the hyperbolic regime as £ crosses the bifurcation point, (d) The 
velocity along the cylinder's axes displays small oscillations, but becomes more negative under 
the influence of the oscillations in the 6-t; plane. 



5.3. (1, l)^- planar periodic orbits 

5.3.1. Existence Orbits of the plane with unit ratio between the y and z velocity components 
are peculiar in that they remain periodic even when, for h > 0, they hit the inner walls of the 
cylinders, i.e. at y,z = ±h. This is generally not so of (l,n) PPOs, which typically cease to 
exist after a collision with the inner walls of the cylinders. 

We define the (1,1)/; planar periodic orbits to be those among such orbits which make 
exactly k collisions with the inner walls of each cylinder over a periodic cycle. Examples are 



shown in figure 12 
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Figure 12. Families of (1, 1)^ PPOs shown in the plane of the billiard, (a) k = 1, (b) k = 2, (c) 

k = 3. 



We identify these orbits by the initial coordinates 
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with the corresponding velocity 
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The values of the geometrical parameter, /i, for which these orbits can be observed are 
restricted to the interval h m [ n < h < h max , whose bounds correspond respectively to the value 
of h for which the orbit hits the corners at the intersection of outer cylinders walls (in which 
case the segments of the orbits have constant lengths), and to the value of h for which the 
longest segment grazes the corner at the intersection of the inner cylinder walls: 

h min =2(fc-l), (2la) 
h max = 2(k+l). (21/7) 

Given a value of h in this interval, the dynamical parameter, £, which characterises the 
orbits of this family, is itself restricted to the interval £ m i n < £ < £ m ax, where 

£ min =max[/i-2(£+l),2(£-l)-/z], (22a) 
£ max = min[2(fc + 1) - h, h - 2{k - 1)] . {22b) 

Interestingly, the (1, l)k+\ orbit at h — h m [ n is identical to the (1, 1)^-1 orbit at h = h mSLX , 
except for eight segments which connect the inner to the outer corners back to back. 



5.3.2. Stability The stability analysis of the (1 , 1)^ PPOs proceeds along the lines of sections 
5.2.2| and |5.2.3[ The linear map is of the form ( fT2| ) and the quadratic terms similar to equations 
( |17a[ )-( [T7^1 ) and ( |18a[ )-( [T8^| ). The coefficients are, however, different and their expressions 
get more complicated as k increases. 




Figure 13. Results of the combined linear and nonlinear stability analysis of the (1, 1)^ PPOs. 
for k = 1,2,3. The conventions are the same as in figure [9] 
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The results of this analysis are shown in figure [13} The patterns are very similar to 
those displayed in figure [9] for the (l,l)o PPOs, with the conclusion that elliptic regions 
recur for every value of k at geometrical parameter values of h near h = 2k. To be more 
precise, the elliptic regions come in the shape of two tongues, elongated along the lines 
h = 2k — 1/(2/:+ 1) ± £, which are symmetric reflections of each other along the line 
h — 2k— l/(2k + 1). Moreover, the lower tongues are nonlinearly stable around £ = in 
the interval 2k - 2 /(2k + 1) < h< 2k - I /(2k + 1) and the upper tongues unstable. 

5.4. (m,n)o~ planar periodic orbits 

5.4.1. Existence The (m — 1 , n)o planar periodic orbits can be identified by the phase point 



with position coordinates qo identical to that of the (1, l)o orbits, equation ( 14), and velocity 



Po = 7 =L= (0,l,n). (23) 
v 1 +^ 

These orbits exist for geometric parameter values between the bounds 

^min = - 2, (24a) 
/*max = -. (24b) 

n 

Given a value of h in this interval, the dynamical parameter, £, is restricted, for n odd, to 
the interval £ m i n < £ < £ m ax, where 

£min = max(n/i — 2, — 2 — h) , (25a) 

£max = min(2 - nh, h + 2) , (25Z?) 
and, for n even and < 0, to the interval £ mm < £ < £ m ax, where 

£min = -(2 + h), (26a) 

£max = 2 + /z, (26b) 

or, for h > 0, to the two intervals £ mm < |e| < £ m ax, where 

£ min =(2n+\)h, (21a) 
£max = 2 + /L (21b) 

5.4.2. Stability The linear stability analysis of the (1 , n)o PPOs can be carried out along the 



lines of the analysis of the (l,0)o PPOs described in section 5.1.2 The result yields a linear 
map of the form ( [T2| ), with coefficients, given, for the n = 2 case in |Appendix A.2[ equations 
( |1.3a| )-( [T3g| ). The results are summarised in figure [14} 



Regions of the parameter space corresponding to elliptic eigenvalues of the linear form 
( [T2| ) are restricted to the h < — 1 area. They are divided into n strips, each extending over the 
whole interval — 2 — h < £ <2 + h. These strips are all bounded from below by a straight line 
and crossed in the middle by a line of parabolic eigenvalues. Their upper bounds are rather 
complicated functions of the parameters. 

We note that, as n increases, the elliptic areas fill out the whole area — 2 < h < — 1, 
— 2 — h < £ <2 + h, resembling the pattern observed for the ( 1 , 0)o PPOs shown in figure [7} 
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Figure 14. Results of the linear stability analysis of the (l,n)o PPOs plotted in the parameter 
space. The conventions are identical to figure [7] 



This is not surprising, since, as n increases, the (l,n)o PPOs can be viewed as perturbations 
in the plane of the billiard of the the (1 , 0) PPOs. 

The interesting feature of the (l,n)o PPOs is that, as opposed to the (l,0)o PPOs, their 
nonlinear stability analysis is not trivial. As with the linear stability analysis, we proceed 
in a way similar to section |5.2.3| to obtain quadratic corrections to the linear stability map 
( [T3] ), similar to equations ( |17a| )-( [T7^| ) and ( |18^| )-( [T8^| ). We provide in section |Appendix| 



A.2 the explicit coefficients of the quadratic forms associated to the (l,2)o PPOs, see 



equations ( |1.4a[ )-( [L4c| ) and ( |1.5a| )-( [T3c| ). The expressions coefficients of the quadratic forms 
associated to the (l,n)o PPOs, n > 3, will not be given explicitly. 



The results of the nonlinear stability analysis are displayed in figure \T5\ and suggest that 
every elliptic strip may give rise to nonlinearly stable oscillations, at least over some range of 
values of the geometric parameter h. It is also interesting to note that the coefficients of the 
quadratic forms are all identically zero at the n values of h which correspond to the straight- 
line lower bounds of the tongues of elliptic stability of these orbits. 



Examples of such stable oscillations are shown from actual trajectories in figure 16 



including the example of nonlinear oscillations around a (2,5)o PPO. This orbit, as well as 
other (m, n)o PPOs, display features similar to those described above for the (1 , n)o PPOs, and 
will we not go into further details. 
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Figure 15. Results of the nonlinear stability analysis of the (l,n)o PPOs, for n = 2,3,4, 
combined with those of the linear analysis shown in figure [15] The conventions are the same 
as in figure [9] 
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6. Helical orbits 

We now proceed to analyse orbits of helical shape, which wind around each cylinder in the 
squared billiard in turn. 

By symmetry, helical orbits of any pitch may occur along the surfaces of infinitely long, 
straight cylinders. When these orbits cross the intersection between two cylinders, however, 
they typically lose focus and go astray as they hit a corner between the two cylindrical 
surfaces. Nonetheless, there are isolated points where the intersection of the surfaces of 
the two cylinders is smooth, and a helical orbit will thus be transmitted unscathed if it 
passes through these specific points. One may further expect these orbits to generate stable 
oscillations in their immediate vicinity. Indeed, so long as the orbits remain close enough to 
the surface of the billiard, with a short separation between successive collisions, no defocusing 
can take place. We emphasise that this phenomenon is specific to higher-dimensional billiards 
with curved surfaces. The problem of assessing the stability of such orbits is again a fully 
nonlinear one which requires a suitable treatment. 

In the case under consideration, where the cylinders' axes are in the x = plane and 
intersect at right angles, the points of smooth intersections are the points x = ±1, where the 
curvature with respect to the other y or z coordinate vanishes. Thus a helical orbit which 
crosses over from one cylinder to the other atx = ±l,y = ±(l+/z) and z = y 9 or z = —y 9 is a 
smooth periodic orbit of the squared cylindrical stadium. (The conditions must be verified at 
every intersection.) These are rather stringent conditions and impose restrictive symmetries 
on the orbits. 

We characterise these helical periodic orbits in terms of the number of times that they 
whirl around specific axes. Let n\ and n2 denote the number of half periods the helices make 
about the cylinders along the y and z axes, respectively. By symmetry, n\ and n^ must be odd 
integers, and it is easily found that the (^i,^2)-helix exists only dXh — ^/n\ni7i/2 — 1, with 
initial condition in the z = plane at x = 0, y = h if (n\ — l)/2 is even, y = h + 2 if odd, 
and with initial velocity components ±y/n\/{n\ +^2) along the x-axis and ±^^2/ (^1 +^2) 
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along the z-axis. Technically, these orbits are of infinite period and are therefore impossible to 
realise. However, they give rise to an infinite sequence of periodic orbits of finite periods 
which approximate them. We will denote these approximate helices by (^i,^2)y^-helical 
periodic orbits (HPO), meaning that each half period (of which there are 2{n\ +^2) along 
the orbit) is made out of k segments. Examples of such long periodic orbits are shown in 
figure 17; see also related media M 




Figure 17. Periodic orbits approximating exact helical periodic orbits of the squared 
cylindrical stadia. They correspond respectively to: (a) a (1, 1) 101 -periodic orbit 
approximating the corresponding helix near h = 7t/2 — 1; (b) a (1,3) 101 -HPO near h = 
V3x/2- 1; (c) a (1,5) 101 -HPO near h = V5x/2-l; (d) a (3, 3) 101 -HPO near h = 3^/2-1; 
(e) a (3,5)ioi-HPO near h = y/\5n/2-\\ and (f) a (5,5)ioi-HPO near h = 5n/2-l. The 
exact values of h where the actual (^1,^2)^=101 periodic orbits can be found are provided in 
the text. Regardless of the actual period, they occur in an interval of values of h for symmetric 
cases such as (a), (d) and (f), but only at isolated values of h for asymmetric cases such as (b), 
(c) and (e). 



An important point is that, irrespective of their stability, finite-period approximations 
to helical periodic orbits can be moved around in some neighborhood, so long as they pass 
through the intersection between two cylinders away from their boundaries and do not hit 
a sharp angle. Finite-period approximations to a helical periodic orbit thus belong to one- 
parameter continuous families of such orbits, indexed by the dynamical parameter £ which 
we define below. 



II 



http : //www.ulb . ac.be/~tgilbert/research/ 
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We now turn to the properties of these orbits and analyze their stability for specific 
classes. 



6.1. (1, l)^ helical periodic orbits 



6.1.1. Existence Finite-length (1,1)& (k G N, k > 2) helical periodic orbits which 
approximate the (1,1) helix cross the plane intersecting pairs of cylinders perpendicularl)|f[ 
Their velocity at the intersections thus has the form (0, ^Tf)* ^he com pl ete orbits can be 
identified as going through the point with Cartesian coordinate^] 



£odd,(28) 
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The values of the geometrical parameter, h, for which these orbits can be observed are 
restricted to the interval h m [ n <h< h max , whose bounds correspond, respectively, to the value 
of h for which the lengths of the segments of the orbit going through the intersections of the 
cylinders shrink exactly to zero, and to the value of h for which these segments have their 
maximal length, grazing the surface of the cylinders at the intersection, 



^min 
^max 



71 



(*-2)sin— -1, 



71 



(* + 2)sin— -1. 



(30a) 
(30fe) 



Note that, for k large, the width of this interval shrinks to zero like 1/k, with lim^ /i max — h m [ n ^ 
2n/k. The limiting value lim^/^ max = lim^/^ m i n = 7i/2 — 1 is the only value of the geometrical 
parameter for which the (1, 1) helix exists. 

Given a value of h in this interval, the dynamical parameter, £, which characterises the 
orbits of this family, is itself restricted to the interval £ m i n < £ < £ m ax, where 



max 



mm 



71 



71 



1 - h + {k- 2) sin --, 1 +h - (k + 2) sin 
l + /*-(fc-2)sin^,-l-/* + (A; + 2)sin^- 



(31a) 
(31*) 



As with h, these bounds correspond to the values of £ for which the periodic orbit either hits 
a corner at the intersection between two cylinders, or, the other way around, collides with the 
surface of the cylinders in the middle of its longest segment connecting the orbit between two 



transverse cylinders. Examples are displayed in figure 18 



f This property is shared with other symmetrical (n,n)k HPOs. 

+ Although we list only one, there are in total four such similar orbits, which are easily obtained from each other 
by symmetry. This property is shared by all helical periodic orbits of any length. 



Higher-dimensional convex billiards with cylindrical shape 



30 



(a) 




(b) 




(c) 




(d) 




(e) 





(f) 



(g) 




(h) 




(i) 




Figure 18. Families of (1, 1)* HPOs found at h = 0.5. As opposed to the exact helix, which 
is unique, its finite-length approximations are parameterised by a continuous parameter in the 



interval (31a)-(31b). We show them for different values of the parameter taken across the 
interval, for (a) k = 2, (b) k = 3, (c) k = 4, (d) k = 5, (e) k = 6, (f) k = 7, (g) k = 11, (h) k = 23, 
(i)Jfc = 41. 



As opposed to the parameter h which fixes the geometry of the billiard, the parameter £ is 
a dynamical one: it changes as perturbations of the periodic orbit are introduced. Just as was 
the case with the planar periodic orbits, displacements of the helical periodic orbit along the z 
axes as the trajectory returns periodically close to its initial condition ( [28] ) can be interpreted 
as changes in the value of the parameter £; the trajectory comes back closer to another periodic 
orbit than the one it was perturbed away from. In particular the trajectory moves away from 
this family of orbits as soon as the value of £ exceeds either of its two bounds, £ m [ n or £ m ax> 
since it ceases to exist beyond them. 
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6.1.2. Stability 



Birkhoff coordinates. To analyze the stability of the (1, 1)# orbits, we change coordinates to 
the Birkhoff coordinates ( |lla| )-( [Tk/| ), here given by 

yo + l+h \ fcodd, 



arctan : 



x 



= sin ( arctan — — Go ) = 
V u o J 

sin ^ + 8, £odd, 
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k even, 



(32a) 
(32*) 
(32c) 
(32d) 



Linear stability analysis. Introducing perturbations 50, 5z and 5w, the linear stability 
analysis yields the symplectic map 



(33) 
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This form is similar to that of the planar periodic orbits, ( |T2| ), but with extra non-zero 
coefficients m z o and m z ^. Its eigenvalues remain similar: there are two unit eigenvalues 
associated to the motion in the z-w plane, and two nontrivial eigenvalues, r/i and r/2, 
associated to the motion in the 6-% plane. These two eigenvalues can be either hyperbolic, 
when they are real, r/i = i]" 1 , or elliptic, when they have a non-zero imaginary part, r/i = 7]|. 
The matrix elements m a p are generally complicated functions of the two parameters h and £. 

The cases k = 2 and /: = 3 are the simplest and yield results shown in figure [19} exhibiting 
features rather similar to the (1,1)/; PPOs, which are common to all (1, 1)^ HPOs. (Explicit 



values of the coefficients m a p in the linear form (33) for these two cases are provided in 
Appendix"^} ) Namely, there are regions in the e-h parameter plane where r/i and r/2 are 
elliptic, which are in the shape of two symmetric tongues. For k — 2 the symmetry line 
corresponds to h = 3\/2/4 — 1, and to h — 1/3 for k = 3. This value keeps growing for larger 
values of k towards h — 7i/2—\. 

Here, as with the (1,1)/; PPOs, we see that the coefficient m zw in the linear form (33) 
depends only on h - see equations ( |2.1g| )-( [23g| ) in section | Appendix B[ Repeating the 



arguments given in section p22\ we would conclude that all these periodic orbits are linearly 
unstable because of the existence of a bifurcation from elliptic to hyperbolic regimes in the 
range of dynamical parameter values the orbit would visit as it translates further away from 
its initial position under the action of Sw. 



Higher-dimensional convex billiards with cylindrical shape 



32 




■1.0 -0.5 0.0 0.5 1.0 



1.5 r 



1.0 



^ 0.5- 



0.0 



-o.5 y 




-1.0 -0.5 



0.5 1.0 



Figure 19. Results of the linear stability analysis of the (1, 1)2 and (1, 1)3 HPOs. The light 



red areas represent the regions outside the range [£ m in,£max], defined in equation (31a)-(31b). 



The grey areas correspond to regions where the eigenvalues of the map ( [33] ) are hyperbolic, 
and white where they are elliptic. The lines which appear in the elliptic regions are lines of 
parabolic eigenvalues. 



Nonlinear stability analysis. There are, however, again nonlinear effects which may affect 
this scenario, due to terms which are quadratic in the small perturbations. The picture here 
is actually slightly simpler than with the planar periodic orbits. Indeed, since the 0, £ and z 
variables pick up nontrivial perturbations at the order of linear contributions, it is enough to 
perform the quadratic expansion for the w axis alone. This yields the quadratic map 



8w^8w + Q w (80,8Z), 
where Q w (86, 8%) is a quadratic form, given by 



Q W {80,8Z)=(S0 8!; 




(34) 



(35) 



similar to ( |18a| ). The coefficients a W9 b w and c w are here again functions of the two parameters 
h and £. 



Recapitulating the arguments given in section |5.2.3[ the sign of the quadratic form ( [35] ) 
is determined by that of the eigenvalues of its matrix. When the eigenvalues are of the same 
sign, Q w will be positive when they are positive and negative when they are negative. If, on 
the other hand, the eigenvalues are of opposite signs, Q w will oscillate, but, on average, will 
take the sign of the eigenvalue that has the largest magnitude. Stable nonlinear oscillations 
can take place in the parameter regions where the eigenvalue spectrum of Q w changes sign 
so as to confine the small oscillations to an interval of the dynamical parameter where the 



eigenvalues of the linear form ( [33] ) are elliptic. 

Figure [20] shows, for the k = 2 and k = 3 HPO families, the eigenvalue spectrum of the 



matrix of the quadratic form ( [35] ), superimposed on the elliptic regions found from the linear 
stability analysis in figure [19} The coefficients of these quadratic forms are given in | Appendix] 
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I (l,l) 2 HPO 




|Bj The color convention is the same as the one used in section \5?23\ in terms of the sign and 
magnitude of their eigenvalues, which are denoted by X\ and A2. We observe that the elliptic 
regions below and above the lines at h = 3V2/4 — 1 for k = 2 and /i=l/3for/: = 3 display 
clearly different patterns. Whereas the upper elliptic regions are everywhere unstable under 
the quadratic perturbations, the lower elliptic regions contain stable regions around £ = and 
h above /z = — 1 + 1/ y/l for k = 2 and /i=l/6for/: = 3 (these bounds correspond to the lowest 
point of intersection between the elliptic and hyperbolic regions at £ = 0). Though the lower 
bound we read from this analysis is sharp, the upper bound appears perhaps less so. More 
precise results might be obtained from a more detailed analysis of the eigenspectrum of Q w , 
but this will not concern us here. Numerical results, such as shown below, support the claim 
that the upper bound is actually sharp too. The main result here is the existence of a region of 
parameter space where the helical periodic orbits under consideration are nonlinearly stable. 
The stability can be checked numerically by perturbing actual trajectories about these 



periodic orbits, as shown in figure 21 



Similar results, obtained for larger values of k, are shown in figure [22} there resorting 
to a numerical computation of the quadratic form Q w on a grid of points of the parameter 
space. That is to say, we fix numerical values of the two parameters h and £, identify the 
corresponding periodic orbit, and perform a perturbation analysis around that orbit, which 
yields numerical values of the coefficients m a p in equation (33 ) and a W9 b w and c w in equation 
( [35] ). As k gets larger, the regions of stability get thinner, but they persist and can be observed. 
For example, we were able to check the existence of nonlinearly stable periodic orbits for 
the family of (1,1) 101 HPO at /z = 0.5707, which belong to class of orbits displayed in 
figure [TT^ a). The same orbit is, however, found to be unstable at h = 0.5708. 
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e e 

Figure 21. Nonlinearly stable oscillations of the (1,1)2 HPOs measured from actual 
trajectories at h = —0.25 (left column) and h = +0.05 (right column). The oscillations 
measured in the 6-t; (resp. w-z) plane are shown on the second (resp. third) line. Notice 
the order of magnitude of the epsilon oscillations as opposed to those of 8w. 




Figure 22. Numerical results of the nonlinear stability analysis of the (1, 1)& HPOs, for k = 
4,. . . ,9. The conventions are the same as in figure [9] 
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6.2. (n,n)k Helical Periodic Orbits 

6.2.1. Existence We can identify, for arbitrary odd integer n, any finite length (n,n)k HPO 



approximating the (n,n) helices, with the same initial conditions (28)-(29). These are 



basically the same orbits, only the number of collisions along the periodic orbit changes 
from 2k to 2nk and the values of the parameters ( |30a| )-( [30^| ) and ( |31a| )- pT^| ) are accordingly 
changed to 

71 
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hmax = (nk + 2) sin 
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71 
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71 71 

£min = max[— l—h + (nk — 2) sin — - , 1 + h — (nk + 2) sin — ] 



71 



71 . 



£max = min[l + h — (nk — 2) sin — - , — 1 — h + (nk + 2) sin 
Some examples of (3,3)jt and (5,5)^ HPOs are shown in figures [23] and 24, respectively. 



(36a) 
{36b) 

(37a) 
(31b) 



(a) 





(b) 

Figure 23. Families of (3, 3)# HPOs found at h = 3.5, with parameter values across the interval 
gTgHSTg: (a) k = 2,(b)k = 3, (c) /c = 4. 
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6.2.2. Stability The problem of determining the stability of (n,n)k HPOs transposes 
verbatim from section \6.l.2\ The computation gets more challenging as the lengths of the 
orbits increase, but similar results can be obtained. This is supported by the results shown in 



figure 25 , obtained for the (3, 3) 2 HPO. 



4.65685 



4.09117 



3.52548 



2.9598 



2.39411 



1.82843 



• AiA 2 >0, A 1 +A 2 >0 

• A 1 A 2 >0 5 A 1 +A 2 <0 
AjA 2 <0, Aj +A 2 >0 
AjA 2 <0, A| +A 2 <0 




-1. -0.6 -0.2 



0.2 0.6 1. 



Figure 25. Nonlinear stability analysis of the (3, 3) 2 HPO. The conventions are the same as in 
figure [9] 



6.3. (n\,n2)k Helical Periodic Orbits 



6.3.1. Existence. Given odd integers n\ and n^, n\ ^ n^ the (^1,^2)^ HPOs which 
approximate the (ft 1,^2) helical periodic orbit can be identified starting with initial positions 
([28]) and velocity components adapted according to the ratio between the numbers of half 
periods to be completed along the two axes y and z: 
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Because of the asymmetry between the y and z motions, there is here, and for any given 
k, a unique value of the geometric parameter at which these orbits are realised, 
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The corresponding interval of values of the dynamical parameter, £, has the bounds 

2^Jriijn\ sin[n/ (2fc)], &even, 
(1 + y/n 2 /ni) sin[7r/(2£)] 5 fcodd, 
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-2 sin [ft/ (2&)], fceven, 
-(1 + V / n 2 /ni)sin[ft/(2^)], fcodd. 



(406) 



We will focus our attention on three groups of such orbits, namely the (1,3)*, the (1,5)^, 
and the (3,5)* HPOs. Specific examples of families of these classes of (ni,H2)/fc orbits are 



shown in figures 26, 27mndG8l respectively. 



(a) 




(d) 



(b) 




(c) 





(e) 




(f) 




Figure 26. Families of (1 , 3)& HPOs found at the values of the geometric parameter ( 39 ), with 
dynamical parameter values across the interval (40Z? )-(40a ): (a) k = 2, (b) k = 3, (c) k = 4, (d) 
& = 5, (e) k = 6, (f) k = 7. 



6.3.2. 
section 



Stability The nonlinear stability analysis proceeds along the lines described in 
However, since each (n\,n2)k HPO exists only for a singular value of the 



6.1.2 



geometric parameter, we can consider the problem of identifying the elliptic regime of the 

as a function of the dynamical parameter £ alone. 



map ( [33] ) and the quadratic form ([34])- 

This computation can be carried out numerically for moderate values of the overall 



period. The results are displayed in figures [29J [30] and [3TJ respectively, where we plot the 
determinant (=^1^2) and trace (=X\ +A2) of the quadratic form throughout the range of 
dynamical parameter £ values where the eigenvalues of the linear stability analysis are found 
to be elliptic. In all of these cases, we observe that the elliptic regime of dynamical parameter 
values is asymmetric with respect to £ = for k even and symmetric for k odd. In these plots, 
candidate stable HPOs would correspond to values of £ for which both the determinant and 
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(a) 




(b) 




(c) 




Figure 28. Families of (3, 5)& HPOs found at the values of the geometric parameter ( 39 ), with 
dynamical parameter values across the interval (|4Q/7[)-([4Qa|): (a) k = 2, (b) k = 3, (c) k = 4. 



the trace vanish, together with a negative derivative of the trace, i.e. going from positive values 
on the left of this stable point, to negative values on its right. This would ensure the existence 
of an oscillating regime of the quadratic form, which might indeed stabilise the given HPO. 

Considering the stability analysis of the ( 1 , 3 ) ^ HPOS , shown in figure 29 , we observe the 
opposite: namely, the points where both the determinant and the trace of Q w vanish correspond 
to positive values of the derivative of the trace. Therefore they are unstable, since oscillations 
around the given HPO will tend to grow away from it. Similar conclusions are drawn for the 



(3,5)^ HPOs whose stability analysis is shown in figure 30 



The case of the stability analysis of the (1,5)* HPOs, shown in figure |30[ is different. 
There are now two points where both the determinant and trace of Q w vanish. The leftmost 
one is of the same kind as the one observed for the (1,3)^ HPOs. However, the rightmost 
one is certainly a stable point, since it appears to correspond to a minimum of the determinant 
of Q w . It must, however, be noted that it arises near the border of the bifurcation point of £ 
which separates the elliptic region from the hyperbolic one. Stable oscillations around this 
point must therefore be confined to a very small region of phase space. We have not been able 
to obtain conclusive evidence of their existence. 
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Figure 29. Nonlinear stability analysis of the families of (1,3)^ HPOs shown in figure 26 
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Figure 30. Nonlinear stability analysis of the families of (1,5)* HPOs shown in figure 27 



7. Summary and perspectives 

The flat surface component of three-dimensional cylindrical billiards complexifies the stability 
analysis of their periodic orbits. In the presence of two elliptic eigenvalues which generate 
linear oscillations in the plane transverse to the cylinder axis, nonlinear effects along the 
neutral directions may act as a restoring force, allowing for stable oscillations in all phase 
space directions. 

Summarizing our results, the three-dimensional cylindrical stadium billiard with a single 
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oblique plane at angle 71 / 4 with respect to its axis displays a large range of different dynamical 
regimes, stable and unstable, which can be analysed in terms of the spectral properties of 
families of planar and helical periodic orbits of this billiard. We classify these regimes in 
terms of the geometry of the expanded square stadium billiard, which is parameterised by its 
height: 

-2 < h < — I Stable oscillations can be observed in large regions of phase space when 



the height is small enough that opposite circular arcs are closer than their diameter. 
These oscillations take place near a number of families of stable planar periodic orbits. 
Nonlinear effects act on these oscillations as a restoring force which for some of these 
families may constrain the oscillations around the periodic orbit located at the center 
of the family. Though most of the periodic orbits are unstable in some regions of the 
parameter space in this range, others remain stable so that there are no parameter values 
in this region which correspond to a fully chaotic regime. 

I <h < — l + l/y/2 All but the ( 1 , 1 )o planar periodic orbit are unstable in this parameter 



range. Though the (l,l)o PPOs still display regions of elliptic regimes, they are 
destabilized by nonlinear effects. Though the (1,1)2 HPOs exist in this range, they 
are unstable. This range therefore corresponds to a fully chaotic regime: the system 
is ergodic with all four Lyapunov exponents non-zero. 



1 + 1 /y/2 <h< 3^2/4 - 1 The first familiy of helical periodic orbits, the (1, 1) 2 HPOs, 
are nonlinearly stable in this range, around the central orbit at £ = 0. This gives rise to a 
small island of stability in an otherwise chaotic phase space. 



3^/2/4—1 <h<l/6 The (1,1)2 HPOs are nonlinearly unstable in this range and the 



(1, 1)3 HPOs are unstable. This is another range of parameter values for which the system 
is ergodic with all four Lyapunov exponents non-zero. 



1/6 < h < 1/3 The (1,1)3 HPOs are nonlinearly stable around £ = 0. 



1/3 ^h < 7t/2— 1 The sequence of (1,1)^ HPOs, k > 4, follow each other closely in this 
range of parameter values, displaying small regions of nonlinear stability around £ = 0. 
They culminate in the (1,1) helix dXh — 7i/2— 1. 



71/ 2—1 <h < 4/3 The system is ergodic with all four Lyapunov exponents non-zero in this 



range, which separates the last of the (1, 1)# HPOs from the region of nonlinear stability 
of the (l,l)i PPOs. 



4/3 < h < 5/3 The (1, l)i PPOs are nonlinearly stable around £ = in this range. 
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A similar succession of fully chaotic and ergodic regimes, and mixed, non-ergodic, 
regimes where the phase space is separated into a single chaotic region and small elliptic 
regions, continues as we keep increasing the height parameter h. Thus the next regions 
of nonlinear stability correspond to the (3,3)^ HPOs which start for k = 2 near h = 3 and 
culminate at h = 3^/2 - 1, followed, in the interval 18/5 < h< 19/5, by the (1, 1) 2 PPOs, 
and so on. Let us also mention the possibility that stable regimes of oscillations exist for 
isolated values of h, such as with the (1,5)^ HPOs. However, we have not been able to 
confirm their existence numerically. 

We note that the stable oscillation regimes described above are not typical of all billiards 
in the larger class of convex, cylindrical stadium billiards, as they can be easily destroyed 
by a change in the angle of the oblique plane or by the insertion of another oblique plane. 
In particular, these stable oscillations do not exist in the billiards obtained from interacting- 
particle models confined by hard walls. The three-dimensional cylindrical stadium billiard 
with two perpendicular planes at angle n/A with respect its axis is thus ergodic and fully 
chaotic in its whole range of parameter values. It is indeed easy to see that there are no stable 
planar periodic orbits and that helical periodic orbits do not exist. A periodic orbit of this 
billiard always has at least one segment long enough to induce defocusing. 

Cylindrical stadium billiards are easily extended to higher-dimensional billiards and 
provide a fertile ground for exploration of regular and chaotic phenomena in higher- 
dimensional cavities. 
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Appendix A. Stability analysis of the Planar Periodic Orbits 

In this appendix we give the expressions for the coefficients of some of the linear and quadratic 
forms obtained in the stability analysis of planar periodic orbits. 



Appendix A.l. (1,1)q PPOs 



We provide in this section the expressions of the coefficients of the quadratic forms ( |18a| )- 
( lSb) for the (1, l)o PPOs used to obtained the results shown in figure |9j 



a w = -l6V2h 6 + 32V2h 5 £ - 80V2h 5 + l6V2h 4 £ 2 + \AAVlh A £ - UAVlh 4 
-6AV2h 3 £ 3 + 3lV2h 3 £ 2 + 22AV2h 3 £ - 1 llV2h 3 + 16 V2/iV - 160y/2h 2 e 3 
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+l36V2h 2 e - 36V2h 2 + 32V2M 5 + 48V2M 4 - 96\ f lhe i - 16V2M 2 

+2SV2he - 4\[2h - 16^2e 6 + 16V2e 5 + 16^2e 4 - 8 V2e 3 

-4V2e 2 + 2V2e, (1.1a) 

b w = l6V2h 6 - 32V2h 5 £ + 96V2h 5 - \6\f2h 4 £ 2 - \mV2h 4 £ + 2\6\/2h 4 
+64V2h 3 £ 3 - 64V2h 3 £ 2 - 288 V2h 3 £ + 224V2/* 3 - l6V2h 2 £ 4 + \92\[2h 2 £ 3 
-mV2h 2 £ 2 - 224\flh 2 £ + \06V2h 2 - 32^2h£ 5 - 32\[2h£ 4 + 160\/2/ie 3 
-32\f2h£ 2 - 12\[2h£ + 20^2h + 16V2~£ 6 - 32^2£ 5 - 8 V2~£ 4 + 32\[2~£ 3 
+2V2e 2 -8V2e + V2, (\.\b) 

c w = - \6\f2h 6 + 32V2h 5 £ - 1 l2V2h 5 + \6\/2h 4 £ 2 + 116^h 4 £ - 3§4\[2h 4 
-64 V2h 3 £ 3 + 96V2h 3 £ 2 + 352\f2h 3 £ - 400V26 3 + 16 V2h 2 £ 4 - 224 V2h 2 £ 3 
+192V2h 2 £ 2 + 312V2h 2 £ - 260V2h 2 + 32V2h£ 5 + 16^2/ie 4 
-224 V2h£ 3 + \A\\f2h£ 2 + \24\f2h£ - 76 V2h - 16 V2e 6 + 48 V2e 5 
-16^2e 4 - 56V2e 3 + 28V2e 2 + 22^2e - 8^2, (1.1c) 



a z = -32h s - 288/? 7 + 128/? 6 £ 2 + 96h 6 £ - 1040/* 6 + 864/*V + 608/? 5 £ 
- 1904/* 5 - 192h 4 £ 4 - 288/iV + 2096/iV + 1440/? 4 £ - 1828/? 4 - 864/? 3 £ 4 

- 1216/* 3 £ 3 + 2144/iV + 1536/* 3 £ - 852h 3 + 128/* 2 £ 6 + 288/iV - 1072/iV 

- 1600/* 2 £ 3 + 144h 2 £ 2 + 696h 2 £ - I58h 2 + 288/ie 6 + 608/ie 5 - 240/ie 4 
-640/ze 3 - 44/ze 2 + 120h£ - lOh - 32e 8 - 96e 7 + 16e 6 + 160e 5 + 60e 4 
-24e 3 -34e 2 + 16e + 4, (1.2a) 

b z = 32h s + 320h 7 - 128/? 6 £ 2 - 64h 6 £ + 1312/? 6 - 960/? 5 £ 2 - 448/* 5 £ + 2816/x 5 
+ 192h 4 £ 4 + \92h 4 £ 3 - 2152h 4 £ 2 - \2\6h 4 £ + 3336/i 4 + 960/iV + S96h 3 £ 3 
-31\2h 3 £ 2 - 1600ft 3 £ + 2096ft 3 - 128ft 2 £ 6 - \92h 2 £ 5 + 1568ft 2 £ 4 
+ 1408ft 2 £ 3 - 2288ft 2 £ 2 - 1048ft 2 £ + 608ft 2 - 320ft£ 6 - 448ft£ 5 
+896ft£ 4 + 832/ze 3 - 464ft£ 2 - 328ft£ + 64h + 32£ 8 

+64£ 7 - 128£ 6 - 192£ 5 + 104£ 4 + 152e 3 +40£ 2 -48£, (1. 2b) 

c z = -32ft 8 - 352k 1 + 128ft 6 £ 2 + 32ft 6 £ - 1616// 6 + 1056ft 5 £ 2 + 224ft 5 £ 
-3984ft 5 - \92h 4 £ 4 - 96h 4 £ 3 + 3440/zV + 608/z 4 £ - 5636/z 4 - 1056/i 3 £ 4 
-448/z 3 £ 3 + 5536/i 3 £ 2 + S32h 3 £ - 4508/z 3 + 1286 2 £ 6 + 96h 2 £ 5 - 2032/i 2 £ 4 
-704h 2 £ 3 +4456h 2 £ 2 + 6S0h 2 £ - 18626 2 + 352h£ 6 
+224h£ 5 - 1552/ie 4 - 44Sh£ 3 + 1564/i£ 2 + 408/i£ - 326/i - 32£ 8 - 32£ 7 
+208£ 6 + 96£ 5 - 356£ 4 - 136e 3 + 150£ 2 + 144£ - 20. (1.2c) 
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The non-trivial coefficients of the linear form ( 12) for the (l,2)o PPO are found to be 

m ee = 32h 6 + 160/r 5 - 643 4 £ 2 - 64h 4 £ + 232h 4 - 256h 3 e 2 - 2563 3 £ 
+323 2 £ 4 + 64h 2 £ 3 - 272/rV - 2723 2 £ - \92h 2 + 963fi 4 
+ 1923£ 3 + 32h£ 2 + 32h£ - 403 + 72e 4 + 144e 3 + 120e 2 
+ 120e + 49, 



-32h 6 - 2243 5 + 643 4 £ 2 - 5843 4 + 3203 3 £ 2 - 6483 3 - 32h 2 £ 4 



+560h 2 £ 2 



1643 2 - 96h£ 4 + 392h£ 2 + 212/?- 72e 4 + 84e 2 + 120, 



m^ e = -323 6 - 96h 5 +64h 4 £ 2 - 83 4 + 192/rV + I36h 3 - 32h 2 £ A 
+H2h 2 £ 2 + 4h 2 - 96h£ 4 - 12h£ 2 - 68/i - 72e 4 - 36e 2 + 20, 



(1.3a) 



(1.33) 



(1.3c) 



m 



32h 6 + 160h 5 - 64h 4 £ 2 + 64h 4 £ + 232h 4 - 2563 3 £ 2 + 2563 3 £ 



+32h 2 £ 4 - 64h 2 £ 3 - 2723 2 £ 2 + 212h 2 £ - \92h 2 + 963£ 4 - \92M 3 
+32h£ 2 - 32h£ - 40h + 72e 4 - 144e 3 + 120e 2 - 120e + 49, 



m 



zw 



20^53 + 40^5, 



(1.3d) 
(\.3e) 



The coefficients of the quadratic forms ( |18a[ )-( T83"| ) obtained for the (l,2)o PPOs are 



4(23 + 3) 



1283 9 - 643 8 (4e - 1 1) - 323 7 (4e 2 + 48e - 33) 



V5 

+ 163 6 (32e 3 - 20e 2 - 182e - 27) 
- 163 5 (8e 4 - 160e 3 - 42e 2 + 46e + 1 15) 

-83 4 (32e 5 + 120e 4 - 472e 3 - 298e 2 - 395e + 10 
+83 3 (l6e 6 - 128e 5 - 292e 4 + 16e 3 + 1 12e 2 + 233e + 159 
+4h 2 (l44£ 6 - 280e 5 - 468e 4 - 780e 3 - 452e 2 - 353e - 1 1) 
+43 (216e 6 + 24e 5 + 108e 4 - 186e 3 - 174e 2 - 171e - 95) 
+432e 6 + 504e 5 + 864e 4 + 828e 3 + 540e 2 + 347e + 100 



(1.4a) 



3 W = 



4(23 + 3) 



1283 9 - 643 8 (4e - 15) - 323 7 (4e 2 + 56e - 81) 



V5 

+ 163 6 (32e 3 - 52e 2 - 286e + 155) 
- 163 5 (8e 4 - 192e 3 + 1 18e 2 + 282e + 69) 

-83 4 (32e 5 + 88e 4 - 824e 3 + 174e 2 - 47e + 405 
+43 3 (32e 6 - 320e 5 - 328e 4 + 1344e 3 + 192e 2 + 848e - 147) 
+23 2 (288e 6 - 1 136e 5 - 376e 4 - 56e 3 + 656e 2 + 600e + 685) 
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+2h (432e 6 - 816e 5 + 168e 4 - 1040e 3 + 26e 2 - 300e + 149^ 
+432e 6 - 360e 5 + 360e 4 - 600e 3 - 282e 2 - 250e - 245 



(1.4b) 



128ft 9 - 64ft 8 (4e - 19) - 32ft 7 (4e 2 + 64e - 145) 



4(2/i + 3) 
Cw ~ V5~ 
+16h 6 (32e 3 - 84e 2 - 390e + 545) 

-16h 5 (8e 4 -224e 3 + 342e 2 + 518e-449) 

-8ft 4 (32e 5 + 56e 4 - 1 176e 3 + 1350e 2 + 301e + 104) 

+8ft 3 (l6e 6 - 192e 5 + 28e 4 + 1328e 3 - 1232e 2 + 615e - 689) 

+ 4h 2 (l44e 6 - 856e 5 + 668e 4 + 724e 3 - 396e 2 + 949e - 61 1) 

+4h (216e 6 - 840e 5 + 924e 4 - 854e 3 + 866e 2 - 137e + 205^ 

+432e 6 - 1224e 5 + 1584e 4 - 2028e 3 + 1740e 2 - 857e + 60o] (1.4c) 

a z = - 1024ft 12 - 10752ft 1 1 + 4096ft 10 £ 2 + 4608ft 10 e - 44800ft 10 + 38400ftV 
+42496ft 9 £ - 88000ft 9 - 6144ft 8 e 4 - 13824ft 8 £ 3 + 135424ft 8 e 2 
+152832ft 8 £ - 58464ft 8 - 50688ft 7 £ 4 - 1 12640ftV + 1941 12h 7 e 2 
+2501 12h 7 e + 67728ft 7 + 4096ft 6 £ 6 + 13824ft 6 £ 5 - 150272ft 6 £ 4 
-352000ft 6 £ 3 - 1 1 104ft 6 £ 2 + 1 14592ft 6 £ + 125576ft 6 + 29184ft 5 £ 6 
+97792ftV - 153920ft 5 £ 4 -486656ft 5 £ 3 -338288ft 5 £ 2 - 184800ft 5 e 
+ 15168ft 5 - 1024ftV -4608ft 4 e 7 + 73472ft 4 e 6 + 261376ft 4 e 5 
+108704ftV - 165728ft 4 e 3 - 268248ftV - 223648ft 4 £ - 74900ft 4 
-6144ft 3 e 8 - 27648ft 3 e 7 + 61632ft 3 e 6 + 304384ft 3 e 5 + 3448 16ft 3 e 4 
+278144ft 3 e 3 + 103776ft 3 e 2 + 1 1 664ft 3 £ - 25470ft 3 - 13824ft 2 e 8 
-62208ftV -40608ft 2 £ 6 + 91872ft 2 £ 5 + 176824ft 2 e 4 + 240128ft 2 e 3 
+ 171824ft 2 e 2 + 97624ft 2 e + 22291ft 2 - 13824fte 8 -62208fte 7 
-98064fte 6 - 95904fte 5 - 78432fte 4 - 22224fte 3 + 1262fte 2 + 10760fte 
+6568ft - 5184e 8 - 23328e 7 - 447 12e 6 - 61344e 5 - 67788e 4 - 52200e 3 
-33009e 2 -16080e-3500 (1.5a) 

b z = 1024ft 12 + 12800ft 11 - 4096ft 10 £ 2 - 2560ft 10 e + 67328ft 10 - 44544ftV 
-28160ft 9 e + 189120ft 9 + 6144ft 8 e 4 + 7680ft 8 e 3 -200960ft 8 e 2 
- 129280ft 8 £ + 287456ft 8 + 56832ft 7 e 4 + 71680ftV - 474944ft 7 e 2 
-313600ft 7 e + 172720ft 7 - 4096ft 6 e 6 - 7680ft 6 e 5 + 213760ft 6 e 4 
+275200ft 6 e 3 - 588064ft 6 e 2 - 405280ft 6 e - 125432ft 6 - 31232ft 5 e 6 
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-58880/zV + 410176/i 5 £ 4 + 549120/i 5 £ 3 - 2544166 5 £ 2 - 198880/i 5 £ 
-258240ft 5 + 1024/i 4 e 8 + 2560/zV - 93952/i 4 e 6 - 180480fcV 
+396512/* 4 £ 4 + 573920/i 4 £ 3 + 226136ft 4 £ 2 + 135120/* 4 £ - 87280ft 4 
+6144ft 3 £ 8 + 15360ft 3 £ 7 - 138176ft 3 £ 6 - 277760ft 3 £ 5 + 120528ft 3 £ 4 
+228800ft 3 £ 3 + 305072ft 3 £ 2 + 213760ft 3 £ + 80300ft 3 + 13824ft 2 £ 8 
+34560ft 2 £ 7 - 96480ft 2 £ 6 - 217440ft 2 £ 5 - 96040ft 2 £ 4 - 102000ft 2 £ 3 
+70200ft 2 £ 2 + 55860ft 2 £ + 56074ft 2 + 13824/ze 8 + 34560ft£ 7 - 21 168ft£ 6 
-73440ft£ 5 - 83568ft£ 4 - 128880ft£ 3 - 49004ft£ 2 - 37620ft£ - 6120ft 
+5184£ 8 + 12960£ 7 +4104£ 6 -4320£ 5 - 16344£ 4 - 34020£ 3 - 20910£ 2 
-18360£-8600 (1.5b) 

c z = -1024ft 12 - 14848ft 11 + 4096ft 10 £ 2 + 512ft 10 £ - 93952ft 10 + 50688ft 9 £ 2 
+5632ft 9 £ - 337344ft 9 - 6144ft 8 £ 4 - 1536ft 8 £ 3 + 270592ft 8 £ 2 + 25856ft 8 £ 
-742752ft 8 - 62976ft 7 £ 4 - 14336ft 7 £ 3 + 809024ft 7 £ 2 + 62720ft 7 £ 
-989424ft 7 + 4096ft 6 £ 6 + 1536ft 6 £ 5 - 273152ft 6 £ 4 - 55040ft 6 £ 3 
+1461664ft 6 £ 2 + 81056ft 6 £ - 664408ft 6 + 33280ft 5 £ 6 + 1 1776ft 5 £ 5 
-648000ft 5 £ 4 - 109824ft 5 £ 3 + 1 573648ft 5 £ 2 + 40096ft 5 £ + 64544ft 5 
- 1024ft 4 £ 8 - 5 12ft 4 £ 7 + 1 10336ft 4 £ 6 + 36096ft 4 £ 5 - 902752ft 4 £ 4 
- 1 14784ft 4 £ 3 + 847784ft 4 £ 2 - 24480ft 4 £ + 478612ft 4 - 6144ft 3 £ 8 - 3072ft 3 £ 7 
+ 190144ft 3 £ 6 + 55552ft 3 £ 5 -727952ft 3 £ 4 -46080ft 3 £ 3 -38176ft 3 £ 2 
- 34768ft 3 £ + 302514ft 3 - 13824ft 2 £ 8 - 6912ft 2 £ 7 + 178272ft 2 £ 6 + 43488ft 2 £ 5 
-299304ft 2 £ 4 + 18816ft 2 £ 3 - 335344ft 2 £ 2 + 1 176ft 2 £ + 719ft 2 - 13824/i£ 8 
-6912ft£ 7 + 85104ft£ 6 + 14688ft£ 5 -32832ft£ 4 + 23184ft£ 3 - 176722ft£ 2 
+16920ft£ - 69272ft - 5184£ 8 - 2592£ 7 + 15768£ 6 + 864£ 5 +9612£ 4 
+5400£ 3 - 30309£ 2 + 6480£ - 21020 (1.5c) 

Appendix B. Stability analysis of the Helical Periodic Orbits 

We provide below the expressions of the non-trivial coefficients of the linear stability matrix 



(33) and quadratic form (35) for the (1,1)2 and (1,1)3 HPOs used to obtained the results 



shown in figures [T9| and [20 



Appendix B.l. ( 1 , 1 ) 2 HPOs 



m ee = 256ft 4 - 704V2ft 3 + 1024ft 3 - 512ft 2 £ 2 + 128 V2ft 2 £ - 21 12\/2ft 2 
+2928ft 2 + IQAVlhe 2 - 1024ft£ 2 + 256 V2he - 384ft£ - 2692 V2h 
+3808ft + 256£ 4 - 128 v^fi 3 + 704 V2£ 2 - 944£ 2 + 264 v^£ - 384£ 
-1284^2+1817, (2.1a) 
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m H = -768/? 4 + 1984 V2h 3 - 3072/? 3 + 1536/? 2 £ 2 - 256 V2h 2 e + 5952 Vlh 2 
-8304/? 2 - \9MV2he 1 + 3072/?£ 2 - 512 V2he + 768M + 7412^2/? 
- 10464/1 - 768e 4 + 256 V2e 3 - 1984 V2e 2 + 2736e 2 - 528 V2e + 768e 
+3444^2-4872, {2.1b) 

m^ e = -256h 4 + 832 V2h 3 - 1024/? 3 + 512/? 2 £ 2 + 2496 V2/? 2 - 3536/? 2 
-832V2/ie 2 + 1024/ie 2 + 3548^2/? - 5024/? - 256e 4 - 832 V2e 2 
+ 1168e 2 + 1884^2-2664 (2.1c) 

= 768/? 4 - 2368 V2/? 3 + 3072/? 3 - 1536/? 2 £ 2 - 128 V2/? 2 £ - 7104 V2h 2 
+ 10000/? 2 + 2368^2/ie 2 - 3072/?£ 2 - 256V2/?£ + 384/?£ - 9788^2/? 
+ 13856/1 + 768e 4 + 128 V2e 3 + 2368 V2e 2 - 3280e 2 - 264V2e + 384e 
-5052^2 + 7145, (2.1J) 

m z0 = -128 V2/? 4 -512V2h 3 + 704/? 3 + 256V2/? 2 £ 2 - 128/? 2 £ - 1464V2/? 2 
+21 \2h 2 + 512 V2/?e 2 - 704/?£ 2 + 192 V2he - 256he - 1904^2/? 
+2692/? - 128 V2e 4 + 128e 3 + 472^2e 2 - 704e 2 + 192 V2e - 264e 
-908^2+1284, (2.1e) 



m z< 



I = 384 V2/? 4 + 1536V2/? 3 - 1984/? 3 - 768 V2/? 2 £ 2 + 256/? 2 £ 
+4152 V2h 2 - 5952/? 2 - 1536^2/ie 2 + 1984/?£ 2 - 3UV2he + 512/?£ 
+5232^2/? - 7412/? + 384 V2e 4 - 256e 3 - 1368 V2e 2 + 1984e 2 
-384V2e + 528e + 2436^2 -3444, (2.1/) 

m zw = 8V2/i + 8V2+16. (2.1g) 



a w = 8192 V2h 6 + 1 63 84 V2h 5 e + 49152V2/? 5 - 77824Z? 5 - 8192V2/? 4 £ 2 
+81920V2/? 4 £ - 126976/? 4 £ + 275456V2/? 4 - 389120/? 4 - 32768 V2h 3 e 3 
-32768 V2/i 3 e 2 + 57344/? 3 £ 2 + 358400 V2/? 3 £ - 507904/? 3 £ + 774144 v^/? 3 
-1094144/? 3 - 8192V2/? 2 £ 4 - 98304 V2/? 2 £ 3 + 155648/? 2 £ 3 
- 122880 V2h 2 £ 2 + 172032/? 2 £ 2 + 747520 V2h 2 e - 10565 12h 2 e 
+ 1220384 V2h 2 - 1725952/? 2 + 16384 V2he 5 - 16384V2/?£ 4 + 20480/?£ 4 
-219136V2/?£ 3 +311296/?£ 3 - 1 80224 V2/?£ 2 + 254464/?£ 2 + 775904\/2/?£ 
- 1097216/?£ + 1023552^2/? - 1447536/? + 8192 V2e 6 + 16384 V2e 5 
-28672e 5 - 13312V2e 4 + 20480e 4 - 15360072e 3 + 216832e 3 
-98784 V2e 2 + 139776e 2 + 321248 V2e - 454328e 

+356688^2-504432, (2.2a) 
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b w = -24576^26 6 - 49152 Vlh 5 e - U1456Vlh 5 + 2252806 5 
+24576 Vlh 4 e 2 - 245760V26 4 e + 311136b 4 £ - 794624 Vlh 4 + 1 1264006 4 
+98304 V2h 3 e 3 + 98304 V2h 3 e 2 - 1556486 3 e 2 - 1050624 V2h 3 e 
+ 14909446 3 e -2195456V26 3 + 31032326 3 + 24576V26 2 e 4 
+294912 Vlh 2 e 3 - 4505606 2 e 3 + 329728 Vlh 2 e 2 - 4669446 2 e 2 
-2168832 V26 2 e + 30656006 2 e - 3397056 Vlh 2 + 48040966 2 
-49152 V26e 5 + 49152 V26e 4 - 696326e 4 + 632832 V26e 3 - 901 1206e 3 
+462848 Vlhe 2 - 6538246c 2 - 2227040 V26e + 31493126c - 2796416^26 
+3954768/? - 24576^2e 6 - 49152 Vie 5 + 77824e 5 + 47104 V2e 4 - 69632e 4 
+436224 V2e 3 - 616192e 3 + 2421 12^2e 2 - 342528e 2 - 912224 Vie 
+ 1290104e -956562^2 + 1352784, (1.1b) 

c w = 73728 Vlh 6 + 147456 Vlh 5 e + 442368 Vlh 5 - 6512646 5 
-73728 Vlh 4 e 2 + 737280V26 4 e - 10936326 4 £ + 2292736 Vlh 4 
-32563206 4 - 294912 Vlh 3 e 3 - 294912 V26 3 e 2 + 411191b 3 e 2 
+3078144 V26 3 £ - 43745286 3 £ + 6221824V26 3 - 87966726 3 
-73728V26 2 £ 4 - 884736^26 2 £ 3 + 13025286 2 £ 3 - 880640V26 2 £ 2 
+ 12533766 2 £ 2 + 6285312 Vlh 2 e - 88849926 2 £ + 9450720 Vlh 2 
-133647366 2 + 147456 V26£ 5 - 147456 V26£ 4 + 2334726£ 4 
- 1 824768 V26£ 3 + 26050566£ 3 - 1171456V26£ 2 + 16573446£ 2 
+6379360V26£ - 90209286£ + 7637440^26 - 108010086 
+73728 Vie 6 + 147456V2e 5 - 208896e 5 - 158720V2e 4 + 233472e 4 
_ \134944\fle 3 + 1745664e 3 - 580896 V2e 2 + 821760e 2 + 2582368 Vie 
-3652056c + 2565300^2 - 3627888. (2.2c) 

Appendix B.2. ( 1 , 1 ) 3 HPOs 

m ee = 103686 4 - 13824/z 3 - 207366 2 £ 2 - 17286 2 £ + 6614h 2 + 138246£ 2 
+H51he - 13446 + 10368e 4 + 1728e 3 - 2016e 2 - 168e + 97 (2.3a) 

m e ^ = -31 1046 4 + 362886 3 + 622086 2 £ 2 - 151206 2 - 362886£ 2 + 26646 
-31104e 4 + 4752e 2 -168 (2.36) 

m^ e = -34566 4 + 51846 3 + 69126 2 e 2 - 28326 2 - 51846e 2 + 6646 
-3456e 4 + 912e 2 -56 (2.3c) 



mtf = 103686 4 - 138246 3 - 207366 2 £ 2 + 17286 2 e + 66246 2 + 138246e 2 
- 1 1526e - 13446 + 10368e 4 - 1728e 3 - 2016e 2 + 168e + 97 (1.3d) 
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m z9 = -5lS4V3h 4 + 6912V3h 3 + l036SV3h 2 e 2 + S64V3h 2 e -33l2V3h 2 
-6912^3/ie 2 - 57673/ie + 672^3/i - 5184V3e 4 - 864 V3e 3 
+ 1008V3e 2 + 8473e-48V3 (2.3e) 

m z | = 15552^3/z 4 - 18144^3/z 3 - 31104^3/i 2 e 2 + 7560^3/i 2 + 18144^3/ie 2 
- 1332^3/i + 15552 V3e 4 - 2376^3£ 2 + 84^3 (2.3/) 

m zw = 8^2/2 + 20^2, (2.3g) 

a w = -746496V2/* 6 + 1492992V2/i 5 e + 1617408 V2/* 5 + 746496V2/i 4 e 2 
-261 2736 V27* 4 £ - 1430784 \flh 4 - 2985984V2/i 3 e 3 - 1244160V2/? 3 £ 2 
+1783296V2/z 3 e + 660096V26 3 + 746496V26 2 £ 4 + 32348 16V2/z 2 £ 3 
+746496 v^V - 592704 V26 2 £ - \61\UVlh 2 + 1492992 V2/z£ 5 
-373248 V2/z£ 4 - 1 1 \9144\[2h£ 3 - 190080V2/z£ 2 + 96048 ^2/i£ 
+22008 V2/i - 746496 V2£ 6 - 622080 Vl£ 5 + 20736 V2£ 4 + 122688 V2£ 3 
+ 17136V2£ 2 -6108V2£-1176V2, (2.4a) 

b w = 2239488 Vlh 6 - 441%916\flh 5 £ - 4478976V26 5 - 2239488 V2h 4 £ 2 
+7464960 V2h 4 £ + 3639 168 Vlh 4 + 8957952 V2h 3 £ 3 +2985984 y/2h 3 e 2 
-4852224V2/i 3 £ - 1534464 \flh 3 - 2239488 \flh 2 £ 4 - 8957952V2/, 2 £ 3 
- 1430784 \flh 2 £ 2 + 1534464 \flh 2 £ + 353592 \flh 2 - 4478976v / 2/i£ 5 
+ 1492992 V2/*£ 4 + 2861568 \flh£ 3 + 290304 \flh£ 2 - 235872 V2/*£ 
-42192^2/? + 2239488 V2£ 6 + 1492992 ^2£ 5 - 217728 ^2£ 4 - 290304 V2~£ 3 
-20520 V2£ 2 + 141 12 \fl£ + 2037^2, (2.4Z?) 

c w = -6718464V2/i 6 + 13436928 V2h 5 £ + 12317184V2/, 5 + 671 8464 V2/, 4 £ 2 
-21275136V2/i 4 £ - 9144576V2/* 4 - 26873856V2/i 3 £ 3 - 671 8464 \flh 3 £ 2 
+ 1 3063680 V26 3 £ + 35 14752 V26 3 + 67 1 8464 V26 2 £ 4 + 24634368 \flh 2 £ 3 
+2239488 V26 2 £ 2 - 3872448 V26 2 £ - 737424 V26 2 + 1 3436928 Vlh£ 5 
-5598720 V2h£ 4 - 70917 l2Vlh£ 3 - 279936 Vlh£ 2 + 553392 Vlh£ 
+80136V2/i - 6718464V2£ 6 - 3359232^2£ 5 + 933120^2£ 4 
+637632 V2£ 3 + 9072 V2£ 2 - 30564 V2~£ - 3528^2. (2.4c) 
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